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Abstract 

A  recently  published  experiment  called  "dual  photography"  exploits  Helmholtz 
reciprocity  by  illuminating  a  scene  with  a  pixilated  light  source  and  imaging  other  parts 
of  that  scene  with  a  camera  so  that  light  transport  between  every  pair  of  source-to-camera 
pixels  is  measured.  The  positions  of  the  source  and  camera  are  then  computationally 
interchanged  to  generate  a  "dual  image"  of  the  scene  from  the  viewpoint  of  the  source 
illuminated  from  the  position  of  the  camera.  Although  information  from  parts  of  the 
scene  normally  hidden  from  the  camera  are  made  available,  this  technique  is  rather 
contrived  and  therefore  limited  in  practical  applications  since  it  requires  access  to  the 
path  from  the  source  to  the  scene  for  the  pixilated  illumination. 

By  radiometrically  modeling  the  experiment  described  above  and  expanding  it  to 
the  concept  of  indirect  photography ,  it  has  been  shown  theoretically,  by  simulation  and 
through  experimentation  that  information  in  parts  of  the  scene  not  directly  visible  to 
either  the  camera  or  the  controlling  light  source  can  be  recovered.  To  that  end,  the 
camera  and  light  source  (now  a  laser)  have  been  collocated.  The  laser  is  reflected  from  a 
visible  surface  in  the  scene  onto  hidden  surfaces  in  the  scene  and  the  camera  images 
collect  how  the  light  is  reflected  from  the  hidden  surfaces  back  to  the  visible  surface. 
The  camera  images  are  then  used  to  reconstruct  the  information  from  the  hidden  surfaces 
in  the  scene.  This  document  discusses  the  theory  of  indirect  photography,  describes  the 
simulation  and  experiment  and  used  to  verify  the  theory  and  describes  techniques  used  to 
improve  the  image  quality,  as  measured  by  modified  modulation  transfer  function. 
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RESTORATION  OF  SCENE  INFORMATION  REFLECTED  FROM 
NON-SPECULAR  MEDIA 

I.  Introduction 

A  photographic  technique  known  as  dual  photography,  which  exploits  Helmholtz 
reciprocity,  allows  for  the  position  of  a  digital  camera  and  a  digital  light  source  to  be 
mathematically  interchanged.  This  mathematical  interchange  was  originally  developed  to 
aide  in  the  rendering  of  computer  generated  scenes  and  enables  the  scene  to  be  "viewed" 
from  the  position  of  the  original  light  source  as  though  "illuminated"  from  the  position  of 
the  original  camera.  In  the  original  work  describing  dual  photography,  the  authors'  main 
example  of  their  work  was  to  "show  how  dual  photography  can  be  used  to  capture  and 
relight  scenes."  (Sen,  et  al.,  2005).  Subsequent  work  by  the  authors  which  include 
Compressive  Dual  Photography  concentrate  on  the  creation  of  adaptive  and  non-adaptive 
algorithms  to  more  efficiently  capture  the  large  amounts  of  data  necessary  to  build  the 
light  transport  matrices  required  for  the  technique  to  work.  (Sen  &  Shheil,  2009). 
Because  the  original  goal  of  dual  photography  was  the  rendering  and  relighting  of 
computer  generated  scenes,  no  attempt  was  made  to  recover  details  from  the  scene  not 
directly  visible  to  either  the  camera  or  the  digitized  light  source.  Additionally,  no  work 
has  been  done  describing  the  quality  of  the  dual  image.  Neither  of  these  oversights 
effected  the  exploitation  of  dual  photography  for  the  authors'  original  intended  purposes. 
Nevertheless,  for  applications  outside  the  computer  graphics  community,  the  recovery  of 
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scene  information  not  directly  visible  to  either  the  camera  or  the  light  source  and  a  metric 
of  the  quality  of  the  dual  image  are  of  considerable  interest. 

In  one  of  the  experiments  from  the  original  dual  photography  paper,  a  playing 
card  was  positioned  such  that  the  face  of  the  playing  card  was  not  visible  to  the  camera 
(Sen,  et  al.,  2005).  A  pixilated  light  source  was  placed  with  a  full  view  of  the  face  of  the 
playing  card  and  a  book  was  placed  so  that  when  a  pixel  illuminated  the  playing  card, 
reflections  from  the  card  could  be  imaged  by  the  camera  after  a  intermediary  reflection 
from  the  book  (see  Figure  1). 


Figure  1.  Original  dual  photography  setup  (Sen,  et  al.,  2005)  (reprinted  with 

permission). 

The  pixels  of  the  projector  individually  illuminated  the  playing  card  and  the 
subsequent  reflections  from  the  card  onto  the  book  were  imaged  by  the  camera.  Using  a 
technique  described  in  the  background  section  of  this  document,  the  projector  was 
converted  to  a  "virtual  camera"  and  the  face  of  the  playing  card  was  revealed  to  be  the 
King  of  Hearts. 

While  the  technique  of  dual  photography  is  effective  for  its  original  purpose,  for 
most  applications  outside  the  field  of  computer  generated  graphics,  there  is  no  reason  to 
attempt  dual  photography  as  described  above.  Simply  put,  if  you  are  able  to  get  a 
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pixilated  light  source  in  a  position  to  directly  view  the  object  of  interest,  it  is  much  easier 
to  position  a  camera  in  that  position  and  image  the  object  directly  instead  of  going 
through  the  complicated  and  data  intensive  process  of  creating  a  dual  image.  There  are 
however,  many  applications  where  discretely  viewing  an  object  hidden  from  direct  view 
of  a  camera  is  of  interest.  Extending  the  concept  of  dual  photography  into  one  of  indirect 
photography,  where  neither  the  camera  nor  the  controlling  light  source  has  a  direct  line- 
of-sight  to  the  object  of  interest  would  open  up  countless  new  opportunities  in  the  field  of 
remote  sensing  and  other  fields  of  study. 

This  document  details  the  development  of  the  radiometric  theory  of  indirect 
photography  and  the  experimental  validation  of  that  theory,  in  which  the  image  of  an 
object  was  recovered  without  either  the  camera  or  the  controlling  light  source  having  line 
of  sight  to  the  object  of  interest.  (Figure  2  (b)  is  an  indirect  photograph  of  (a)  produced 
by  a  co-located  digital  camera  and  light  source  neither  of  which  had  direct  line-of-sight  to 
the  object.  Details  will  appear  in  Chapter  IV) 


(a)  (b) 

Figure  2.  Object  of  interest  (a)  and  its  indirect  image  (b)  created  without  either  the 
camera  or  the  light  source  having  line-of-sight  to  the  object 

This  document  begins  with  a  background  section  where  dual  photography  is 
explained  in  detail.  A  brief  review  of  radiometric  principles  and  a  technique  for 
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quantitatively  describing  image  quality  is  also  included.  The  Chapter  III  fully  develops 
the  radiometric  theory  of  dual  photography  and  then  extends  it  to  one  of  indirect 
photography.  The  Chapter  IV  details  the  three  experimental  setups  used  to  validate  the 
theory  of  indirect  photography  and  the  resulting  data.  A  brief  conclusion,  including 
follow-on  research  potentials  is  also  included. 

While  the  technique  of  indirect  photography  is  still  in  the  early  stages  of 
development  and  requires  further  research  before  an  operational  system  can  exist,  the 
ability  to  "see  around  comers"  and  image  hidden  objects  will  have  a  profound  effect  on 
the  intelligence  community.  This  document  lays  the  foundation  for  the  development  of 
that  capability. 
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II.  Background 


This  chapter  reviews  the  concept  of  dual  photography  in  detail  and  provides  a 
description  of  the  algorithm  used  to  produce  a  dual  image.  It  also  reviews  the  basic 
concepts  of  radiometry  and  the  bi-directional  reflection  distribution  function  (BRDF) 
which  are  used  in  the  following  chapters  to  develop  the  concept  of  indirect  photography. 
The  modulation  transfer  function  (MTF)  as  it  relates  to  image  quality  and  concepts  of 
linear  algebra  are  also  discussed. 


Helmholtz  Reciprocity 


Helmholtz  reciprocity  in  its  most  general  form  states  that  the  flow  of 
electromagnetic  radiation,  in  particular  light,  can  be  reversed  without  altering  its  transport 
properties.  Although  in  his  1856  work  on  physiological  optics,  von  Helmholtz  only 
discussed  specular  reflections,  Rayleigh's  1900  work  extended  the  theory  of  reciprocity  to 
include  non-specular  reflections.  In  von  Helmholtz's  own  words  as  quoted  by  Veach 
(Veach,  1998)(von  Helmholtz  &  Southall,  1962): 

Suppose  light  proceeds  by  any  path  whatever  from  point  A  to 
another  point  B  undergoing  any  number  of  reflections  or  refractions  en 
route.  Consider  a  pair  of  rectangular  planes  a  i  and  a2  whose  line  of 
intersection  is  along  the  initial  path  of  the  ray  at  A;  and  another  pair  of 
rectangular  planes  bj  and  b2  intersecting  along  the  path  of  the  ray  when  it 
comes  to  B.  The  components  of  the  vibrations  of  the  aether  particles  in 
these  two  pairs  of  planes  may  be  imagined.  Now  suppose  that  a  certain 
amount  of  light  J  leaving  the  point  A  in  the  given  direction  is  polarized  in 
the  plane  aj  and  that  of  this  light  the  amount  K  arrives  at  the  point  B 
polarized  in  the  plane  b/:  then  it  can  be  proved  that,  when  the  light  returns 
over  the  same  path  and  the  quantity  of  light  J  polarized  in  the  plane  bj 
proceeds  from  the  point,  the  amount  of  this  light  that  arrives  at  the  point  A 
polarized  in  the  plane  aj  will  be  equal  to  K. 
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Apparently  the  above  proposition  is  true  no  matter  what  happens 
to  the  light  in  the  way  of  single  or  double  refraction,  reflection, 
absorption,  ordinary  dispersion  and  diffraction,  provided  that  there  is  no 
change  in  its  refrangibility,  and  provided  it  does  not  traverse  any 
magnetic  medium  that  affects  the  position  of  the  plane  of  polarization,  as 
Faraday  found  to  be  the  case 

Rayleigh's  work  expanded  Helmholtz's  work  to  include  non-specular  reflections.  In  his 
words  (Rayleigh  &  Strutt,  1900): 

Suppose  that  in  any  direction  (i)  and  at  any  distance  (r)  from  a 
small  surface  (S)  reflecting  in  any  manner  there  be  a  situated  a  radiant 
point  (A)  of  given  intensity  of  the  reflected  vibrations  at  any  point  (B) 
situated  in  direction  e  and  at  distance  r'  from  S.  The  theorem  is  to  the 
effect  that  the  intensity  is  the  same  as  it  would  be  at  A  if  the  radiant  point 
were  transferred  to  B.  [Footnote:  I  have  not  thought  it  necessary  to  enter 
into  questions  connected  with  polarization,  but  a  more  particular 
statement  could  easily  be  made.] 

Using  modem  terminology,  if  a  small  reflective  surface  is  illuminated  by  a  small 
light  source,  and  a  small  sensor  is  placed  so  that  it  measures  the  flux  being  reflected  from 
the  surface,  then  the  positions  of  the  light  source  and  the  sensor  can  be  exchanged,  but  the 
measured  reflected  flux  at  the  sensor  will  remain  constant  (Veach,  1998).  It  is  from  this 
concept  of  the  reversibility  of  the  flow  of  electromagnetic  radiation  that  the  concept  of 
dual  photography  was  derived. 


Dual  Photography 

In  a  2005  paper,  Dual  Photography,  Sen  et  al.  "presented  a  novel  photographic 
technique  called  dual  photography,  which  exploits  Helmholtz  reciprocity  to  interchange 
the  lights  and  camera  in  a  scene."  (Sen,  et  al.,  2005)  The  basic  premise  of  dual 
photography  is  to  use  a  pixilated  light  source  to  illuminate  a  scene  of  interest  one  pixel  at 
a  time  and  record  the  reflections  either  directly,  or  after  being  reflected  from  a  second 
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surface.  These  reflections  can  then  be  used  to  create  a  matrix  which  characterizes  the 


light  transport  from  the  light  source  to  the  camera.  This  matrix  can  then  be  transposed 
creating  the  matrix  that  characterizes  the  light  transport  from  the  camera  to  the  light 
source. 

The  authors  of  the  original  dual  photography  paper  used  the  term  "primal 
configuration"  to  describe  the  real-world,  i.e.  physical,  set-up  used  to  record  the  data  and 
the  term  "dual  configuration"  to  describe  the  situation  where  the  camera  and  the  light 
source  have  been  reversed,  in  effect  creating  a  virtual  light  source  out  of  the  original 
camera  and  a  virtual  camera  out  of  the  original  light  source.  These  conventions  will  be 
used  for  the  rest  of  this  document  and  are  explained  in  further  detail  in  the  next  two 
sections. 

Primal  Configuration 

In  the  primal  configuration,  a  projector  with  p  x  q  pixels  is  used  to  light  a  scene 
one  pixel  at  a  time  with  the  resulting  reflections  imaged  by  a  camera  with  a  resolution  of 
mxn,  (see  Figure  3). 


scene 


Figure  3.  Dual  photography  primal  configuration  (Sen,  et  al.,  2005)  (reprinted  with 

permission) 
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Due  to  linearity  of  the  light  transport,  the  system  in  Figure  3  can  be  modeled  as 


(Sen,  et  al.,  2005): 


c'  =  Tp'  (1) 

where  c'  is  an  mn  x  1  column  vector  representing  the  image  recorded  by  the  camera,  p' 
is  an  pq  x  1  column  vector  representing  the  light  source  and  T  is  an  mn  x  pq  matrix 
which  represents  the  light  transport  characteristics  from  each  pixel  in  the  light  source  to 
each  pixel  in  the  camera  (Sen,  et  al.,  2005). 


Dual  Configuration 

Based  on  Helmholtz  reciprocity,  it  is  possible  to  mathematically  interchange  the 
pixilated  light  source  and  the  camera  in  Figure  3  without  altering  the  path  taken  by  the 
light  or  the  energy  transfer  (Sen,  et  al.,  2005).  The  dual  configuration,  shown  in  Figure  4, 
is  one  in  which  the  light  source  and  the  camera  are  mathematically  interchanged  and  can 
be  modeled  as: 

p"  =  Tr  c"  (2) 

where  c"  is  a  mn  x  1  column  vector  representing  the  virtual  light  source,  p"  is  a  pq  x  1 
column  vector  representing  the  virtual  camera,  and  the  matrix  Tt  transpose  of  the  matrix 
T  (Sen,  et  al.,  2005). 
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Figure  4.  Dual  photography  dual  configuration  (Sen,  et  al.,  2005)  (reprinted  with 

permission) 


Dual  Photography  Algorithm 

When  applying  the  concept  of  dual  photography,  for  each  pixel  in  the  light  source 
of  the  primal  configuration,  an  individual  data  image  is  recorded  resulting  in  p  x  q 
images  with  a  resolution  of  mxn.  This  four-dimensional  (4-D)  set  of  data  ( m  x  n  x 
p  x  q)  is  then  "unfolded"  by  creating  a  column  vector  from  each  data  image  and  placing 
it  in  a  column  indexed  to  the  given  pixel  position  to  form  the  2-D  mn  x  pq  transport 
matrix,  T,  which  describes  the  light  transport  characteristics  from  the  light  source  to  the 
object  and  from  the  object  to  the  camera  (see  Figure  5). 

Once  the  transport  matrix  is  known,  Sen  et  al.  intended  for  it  be  used  to  add 
shadows  to  scenes  and  decrease  the  number  of  calculations  required  for  advanced  lighting 
techniques  by  modifying  the  p'  vector  (Sen,  et  al.,  2005)(Sen  &  Shheil,  2009). 

Given  Helmholtz  reciprocity,  the  transpose  of  the  transport  matrix  ,Tt,  can  be 
used  to  describe  the  light  transport  characteristics  from  the  camera  to  the  object  and  from 
the  object  to  the  light  source  (Sen,  et  al.,  2005).  When  Tr  is  multiplied  from  the  right  by 
an  mn  x  1  column  vector,  the  resulting  column  vector  can  then  be  reassembled  to  form 
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the  dual  image  which  has,  in  effect,  turned  the  original  light  source  into  a  virtual  camera 
(see  Figure  6). 


Figure  5.  Creating  the  transport  matrix  from  data  images 


mxn 


Figure  6.  Creating  the  dual  image  from  the  transport  matrix 


The  above  discussion  was  based  on  a  direct  path  from  the  object  to  the  camera. 
This  however,  is  not  a  necessary  condition  for  the  dual  photography  technique  to  be  used. 
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The  authors  of  the  original  paper  demonstrated  this  concept  using  the  set-up  in  Figure  1 
where  the  dual  image  of  the  playing  card  revealed  it  to  be  the  King  of  Hearts  even  though 
the  camera  used  to  create  the  transport  matrix  did  not  have  direct  line-of-sight  to  the 
playing  card. 

Radiometry 

Radiometry  is  the  quantitative  study  of  the  transfer  of  light.  The  following  two 
sections  outline  the  basic  radiometric  principles  necessary  for  understanding  the  non- 
specular  reflections  discussed  later  in  this  document. 

Solid  Angle 

The  basic  unit  of  reflectance  is  the  solid  angle,  which  is  a  3-D  cone  measured  in 
steradians.  The  2-D  analogy  of  the  steradian  is  the  radian  (see  Figure  7).  The  solid  angle 
is  defined  as  (Driggers,  Cox,  &  Edwards,  1999,  p.  91): 

Q  =  ^  (sr)  (3) 

where  A  is  the  surface  area  of  the  sphere  subtended  by  the  solid  angle,  Q  and  r  is  the 
radius  of  the  sphere.  In  much  the  same  way,  the  radian  is  defined  as: 

0  =  ^  (rad)  (4) 

where  s  is  the  length  of  the  arc  subtended  by  the  angle  0  and  r  is  the  radius  of  the  circle. 
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(a) 


(b) 


Figure  7.  Comparing  Radians  and  Steradians 

While  the  area,  A,  is  the  surface  area  of  the  sphere  subtended  by  Q  if  A  «  r 2 
then  A  «  Ad  where  Ad  is  the  area  of  a  flat  plate  that  subtends  the  same  solid  angle  as  A, 
and  Eq.  (3)  can  be  approximated  by  (Dereniak  &  Boreman,  1996,  p.  40): 

Ad 

Q  =  7T  (sr)  (5) 

If  the  normal  to  the  surface,  Ad,  is  not  directed  to  the  vertex  of  the  cone  of  the  solid 
angle,  then  the  projected  area  must  be  used  and  Eq.  (5)  becomes  (Dereniak  &  Boreman, 
1996,  p.  41): 

Ad 

Q  =  —^  cos  9  (sr)  (6) 

where  9  is  the  angle  between  the  surface  normal  and  the  central  ray  from  the  vertex  to 
Ad.  Likewise,  the  differential  solid  angle  dQ  can  be  written  as  (Dereniak  &  Boreman, 
1996,  p.  44): 

dAd 

dQ.  =  — —  cos  9  (sr)  (7) 
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Radiometric  Quantities 


The  basic  radiometric  equation  is  flux  radiated  per  projected  unit  area  of  the 
source  per  solid  angle  of  the  detector,  mathematically  described  by  (Dereniak  & 
Boreman,  1996,  p.  45): 


Le  =  — - 7— (Wcm  2sr  x) 

dAs  cos  9S  d£ld 


where  Le  is  the  radiance  of  the  surface,  d<De  is  the  differential  flux,  dAs  cos  9S  is  the 
differential  projected  area  of  the  source,  and  dQd  is  the  differential  solid  angle  subtended 
by  the  detector.  The  flux  in  radiometric  equations  is  typically  quantified  in  one  of  two 
units:  mks/Joules  denoted  by  a  subscripted  e,  or  photon  units  denoted  by  a  subscripted  p. 
While  either  is  valid,  mks  units  will  be  used  for  this  entire  document. 

All  other  radiometric  quantities  can  be  derived  from  this  basic  equation.  By 
rearranging  Eq.  (8)  to  isolate  flux,  it  becomes  (Dereniak  &  Boreman,  1996,  p.  45): 


cos  6S  dAsdQ.d(W) 


To  obtain  intensity,  which  is  flux  per  solid  angle,  Eq.  (9)  becomes  (Dereniak  &  Boreman, 
1996,  p.  46): 


which  can  also  be  written  in  differential  form: 

dle  =  Le  cos  9S  dAs(\IV  sr-1) 


(10) 


(11) 
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Likewise,  to  obtain  exitance,  which  is  flux  per  unit  area  of  the  source,  Eq.  (9)  now 
becomes  (Dereniak  &  Boreman,  1996,  p.  46): 


Me  =  (^)  =  J  Le  cos  es  dQd(W  cm  2)  (12) 

which  can  also  be  written  in  differential  form: 

dMe  =  Le  cos  9S  dCld(yV  cm-2)  (13) 

The  last  equation  typically  used  in  radiometry  is  also  the  only  one  from  the  point 
of  view  of  the  detector.  Irradiance  is  the  flux  per  unit  area  incident  on  the  detector. 
Using  Eq.  (6),  the  differential  solid  angle,  dQd,  can  be  rewritten: 

dQd  =  cos  Qd  (14) 

where  r  is  the  range  between  the  source  and  the  detector.  Substituting  Eq.  (14)  into 
Eq.  (9)  yields: 

Oe  =  Jj  Lg  cos  es  dAs  cos  (15) 

By  combining  the  r2  with  dAs  cos  ds  into  the  differential  solid  angle  subtended  by  the 
source: 


Eq.  (15)  becomes: 


dAs 

d£ls  =  -pj-  cos  9S 


Ot 


T 

Jj 


LP  cos  9, 


dQsdAs 


(16) 


(17) 


14 


and  the  irradiance  becomes  (Driggers,  Cox,  &  Edwards,  1999,  p.  92): 


Ee 


Le  cos  9d  d£2s(W  cm  2) 


or  in  differential  form: 


(18) 


dEe  =  Le  cos  9d  d0.s(W  cm  2)  (19) 

In  the  equations  above,  if  the  areas  of  the  source  or  detector,  As  or  Ad,  are  small 
compared  to  the  range  squared,  r2,  the  small-angle  approximation  can  be  invoked  for  a 
uniform  source.  It  assumes  the  irradiance  on  the  detector  can  also  be  considered  uniform; 
therefore,  the  radiometric  equations  can  be  simplified  to: 


0e  =  Le  cos  0S  AsQd(W) 

(20) 

4  = 

/d®e\  . 

1  g^J  =  Le^sesAs(yvsr-1) 

(21) 

Me  =  \ 

/dOe\  , 

ldi4  J  =  EecosdsQd(W  cm  2) 

(22) 

Ee  =  [ 

^®e\ 

ydAd)  =  Le  C°S  6d  Qs(^W  Cm  ^ 

(23) 

A  summary  of  the  basic  radiometric  quantities  is  contained  in  Table  3.  These 
equations  will  be  used  to  create  the  bi-directional  distribution  function  in  the  next  section. 
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Table  3.  Basic  Radiometric  Quantities 


Symbol 

Quantity 

Radiometric  Ratio 

Units 

Qe 

Energy 

J 

Flux 

W 

Le 

Radiance 

32®e 

dAs  cos  9S  dQ.d 

d®e 

W  cm"2sr_l 

W  sr'1 

h 

Intensity 

d£id 

d<t>„ 

Me 

Exitance 

fa's 

d<t>e 

W  cm'2 

Ee 

Irradiance 

dA~d 

W  cm'2 

Bi-directional  Reflectance  Distribution  Function 


The  bi-directional  reflectance  distribution  function  (BRDF)  was  initially 
described  by  Nicodemus  in  1977,  and  is  generally  defined  as  the  ratio  of  the  radiance 
reflected  from  a  surface  to  the  irradiance  incident  on  the  same  surface  from  a  given  solid 
angle  (see  Figure  8)  (Nicodeums,  Richmond,  Limperis,  Ginsberg,  &  HSIA,  1977) 


BRDF (0, 0, 9 0  ,  x,  X) 


dLe(9, 0, 0,0,  x, X) 
dEe(9, 0,  x,  X) 


(sr  x) 


(24) 


or  alternately  (Stover,  1995,  p.  21): 


f  (.9, 0, 9 ,  (p ,  x,  X) 


dLe(9, 0, 0,0,  x,  X) 
dEe(9,cj),x,X) 


C sr  -1) 


(25) 


where  9  and  9 '  are  the  respective  incident  and  reflected  elevation  angles  with  respect  to 
the  surface  normal,  0  and  0  are  the  respective  incident  and  reflected  azimuth  angles  with 
respect  to  a  coordinate  system  in-plane  with  the  reflecting  surface,  x  is  the  position  on  the 
reflector  and  A  is  the  wavelength  of  the  radiation  (see  Figure  8).  Polarization  can  be 
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handled  using  a  Stokes  vector  to  represent  the  incident  irradiance,  and  elements  of  the 
Mueller  matrix  are  individual  BRDFs  (Bickel  &  Videen,  1991)(Pezzaniti,  Chipman,  & 
McClain,  1994).  Furthermore,  Eq.  (25)  can  be  rearranged  to  give  the  radiance  from  a 
surface  given  the  irradiance  and  the  BRDF  of  that  surface: 

dLe(0,  <p,  O',  <p',x,X)  =  dEe(9,<p,x,X)f(d,<p,d',<p  ,x,X)(yVcm~2sr~1)  (26) 

which  will  be  used  to  develop  the  radiometric  theory  of  dual  and  indirect  photography. 


Figure  8.  BRDF  reference  angles  (Balling,  2008) 


Glint  Angle 

The  glint  angle  is  a  construct  used  to  model  the  BRDF  of  surfaces.  Simply  put, 
the  glint  vector  bisects  the  incident  irradiance  vector  and  reflected  radiance  vector,  in  the 
plane  formed  by  the  vectors.  The  glint  angle  is  the  elevation  angle  of  the  glint  vector 
with  respect  to  the  surface  normal  (Sundberg,  Gruninger,  Nosek,  Burks,  &  Fontaine, 
1997).  See  Figure  9  where  G  is  the  glint  vector  and  0  is  the  glint  angle.  In  the  model 
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chosen  to  simulate  the  dual/indirect  photography  experiments,  reflected  vectors  which 
produce  the  same  glint  angle  are  assigned  the  same  radiance  (Torrance  &  Sparrow,  1967) 
(Beard  &  Maxwell,  1973). 


Naz 


Figure  9.  Glint  vector  and  corresponding  glint  angle 


Micro-facet  BRDF  Model 

The  micro-facet  model  assumes  the  reflecting  surface  is  comprised  of  a  collection 
of  small  micro-facts  each  obeying  Snell's  law  of  reflection.  Each  micro-facet  is 
characterized  by  its  local  surface  normal  and  the  micro-facets  are  then  distributed, 
generally,  symmetrically  about  the  global  surface  normal.  A  well  studied  form  of  the 
BRDF  using  the  glint  angle  is  given  by  (Priest  &  Meier,  2002): 


(27) 


where  0  is  the  glint  angle  as  described  above  which  can  be  found  using  Eq.  (28)  (Priest  & 


Meier,  2002): 


18 


(28) 


cos(0)  = 


cos  ( 9i )  +  cos  (0r) 
2  cos  (/?) 


where  (Priest  &  Meier,  2002): 


cos(2/?)  =  cos(0j)  cos(0r)  +  sin (6^)  sin(0r)  cos  (A <p)  (29) 

and  is  the  incident  elevation  angle  with  respect  to  the  global  normal,  Qr  is  the  reflected 
elevation  angle  with  respect  to  the  global  normal,  and  A0  is  the  difference  between  the 
incident  and  reflected  azimuthal  angles  (see  Figure  8). 


Modulation  Transfer  Function 


The  quality  of  an  image  can  be  characterized  in  two  ways:  1)  the  impulse 
response  of  the  system  or  2)  its  Fourier  Transform,  the  optical  transfer  function.  The 
impulse  response  which  is  also  known  as  the  point  spread  function  (PSF),  is  the  2-D 
response  of  the  system  to  a  point  source  (Dereniak  &  Boreman,  1996,  p.  505)(Gaskill, 
1978,  pp.  334-336). 

The  object  recorded  by  an  imaging  system  can  be  described  by  its  radiance  as  a 
function  of  position  #(x,y).  This  function  can  be  further  broken  down  into  a  series  of 
evenly  spaced  point  sources  with  intensities  equal  to  the  intensity  of  the  object  at  that 
point  (Goodman,  2005,  p.  21).  Assuming  the  system  is  linear  shift  invariant  (LSI),  the 
PSF  of  each  point  of  the  object,  A(x,  y),  on  the  image  plane,  can  be  summed  to  form  the 
image,  g>(pc,y').  Another  way  of  describing  this  model  of  an  image  is  the  convolution  of 
the  object  with  the  PSF  of  the  imaging  system  (Dereniak  &  Boreman,  1996,  p.  506). 
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&(x,y)  =  fix, y)  **A(x,y) 


(30) 


where  **  denotes  a  2-D  convolution. 

By  using  the  convolution  theorem  for  Fourier  transforms,  Eq.  (30)  can  be 
rewritten  (Goodman,  2005): 

Get, o=  Ftf.o  *m. o  (3i) 

where 

=  (32) 

n(f,0  =  mMx,y» 

and  $  is  the  Fourier  Transform,  <f  is  the  spatial  frequency  in  the  x  direction  and  (  is  the 
spatial  frequency  in  the  y  direction.  The  function  H,  which  is  also  known  as  the  Optical 
Transfer  Function  (OTF),  describes  the  ability  of  the  system  to  transfer  the  object's 
spatial  distribution  of  light  to  the  image  plane  (Dereniak  &  Boreman,  1996,  p.  507).  The 
OTF  is  generally  a  complex  valued  function;  therefore,  it  can  be  described  with  both 
amplitude  and  a  phase  (Dereniak  &  Boreman,  1996,  p.  507)  (Gaskill,  1978,  p.  358). 

OTF O  =  \HiS,0\  exp  [WiS,  0}  (33) 

The  modulus  of  the  OTF,  |H(<f,  <■)!,  is  the  modulation  transfer  function  (MTF) 
and  describes  the  imaging  system's  ability  to  transfer  the  spatial  frequency  of  light  to  the 
image  plane.  Likewise,  the  argument  of  the  OTF,  0(f,  O  is  the  phase  transfer  function 
(PTF)  and  describes  the  phase  response  due  to  an  asymmetry  in  the  PSF  (Wolfe  &  Zissis, 
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1993,  pp.  8-31).  Based  on  Eq.  (31),  one  way  to  find  the  MTF  without  explicitly  knowing 
the  OTF  or  the  PSF  is  (Goodman,  2005,  p.  139): 


G(U) 

Ftf.O 


(34) 


Another  way  to  think  of  the  MTF  is  as  a  measure  of  the  relationship  between  the 
brightest  and  dimmest  portions  of  the  image  with  that  of  the  average  level.  It  is  described 
in  the  IR  Handbook  as  a  measure  of  what  is  commonly  referred  to  as  contrast  (Wolfe  & 
Zissis,  1993,  pp.  8-31): 


Lmax  '  Lmin  (fvh) 

Lmax  (?i.Ci)  +  Lmin 


(35) 


where  Lmax  is  the  maximum  radiance  and  Lmin  is  the  minimum  radiance  at  the  specific 
spatial  frequencies  ^  and  (1.  Eq.  (35)  is  most  often  used  when  the  object  is  sinusoidal  or 
has  regularly- spaced  bars  such  as  those  commonly  found  in  Air  Force  bar  charts.  Due  to 
the  nature  of  the  objects  chosen  for  the  experiment,  Eq.  (35)  will  be  modified  to: 


White  —  Black 
White  +  Black 


(36) 


where  White  is  the  average  radiance  of  all  the  pixels  in  the  test  image  that  are  white  in 
the  ideal  object  and  Black  is  the  average  of  the  pixels  in  the  test  image  that  are  black  in 
the  ideal  object.  The  MTF  described  above  will  be  used  to  quantitatively  describe  the 
ability  to  resolve  spatial  frequencies  in  both  dual  and  indirect  images1. 


1  Note:  While  a  traditional  MTF  ranges  from  0,  (no  modulation),  to  1,  (no  decrease  in  modulation  from  the 
object),  the  modified  MTF  ranges  from  -1,  (a  perfect  negative  of  the  image),  to  0  a  uniform  i.e.  gray  image 
to  1  a  perfect  replication  of  the  image. 
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Linear  Shift  Invariant  Systems 

For  Eq.  (31)  to  hold,  the  system  must  be  a  Linear  Shift  Invariant  (LSI)  system 
(Gaskill,  1978,  p.  139).  A  system  is  said  to  be  linear  if  for  a  system  characterized  by  the 
operator,  L,  then  for  two  arbitrary  signals  /i(x)  and  /2(x)  such  that  (Boas,  1983,  p.  127): 


£{AO)}  =  0i(x) 

£{/2(x)}  =  02  00 


(37) 


and  constants,  ax  and  a2,  then 


£{ai/i(x)  +  a2/2(x)}  =  a1g1  (x)  +  a2g2(x)  (38) 

Likewise,  a  system  is  said  to  be  shift  invariant  if  "the  only  effect  caused  by  a  shift 
in  the  position  of  the  input  is  an  equal  shift  in  the  position  of  the  output"  (Gaskill,  1978, 
p.  139).  Therefore,  a  system  is  shift  invariant  if  given: 

£{/i(*)}  =  0i  00  (39) 

then: 


£{/i(x  -  x0)}  =  g^x  -  x0 ) 

The  PSF  of  the  dual  and  indirect  images  have  been  assumed  to  be  LSI. 


(40) 


Linear  Algebra 

Two  matrix  multiplication  concepts  of  linear  algebra  which  are  used  in  this 
research  are  the  Hadamard  and  Kronecker  products  of  matrix  multiplication.  The 
Hadamard  product,  which  is  sometimes  referred  to  as  entrywise  product,  is  formally 
defined  by  the  following. 
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Given  two  matrices  of  the  same  dimensions. 


A,®  6  Rmxn  (41) 

then  the  Hadamard  product  of  A  and  IB  is  defined  (Hogban,  Brualdi,  Greenbaum,  & 
Mathias,  2003): 

(A  O  1)  ij  =  A itJ  x  Mtj  (42) 

then 


A  ©  E  G  Rmxn 

The  Kronecker  product  is  defined  as  follows: 
Given 


Ae  Rmxn 
IB  G  Rpxq 


(43) 


(44) 

(45) 


then  the  Kronecker  product  of  A  and  IB  is  defined  as  (Hogban,  Brualdi,  Greenbaum,  & 
Mathias,  2003): 


A®  IB 


'a1;L] B  •••  amlIB' 

...  &TJZTI  IB  - 


where 


A  0  1  G  R mpxnq  (46) 

Another  linear  algebra  operator  used  in  this  document  is  the  Vec  operator. 
The  Vec  operator  takes  a  matrix  as  its  input  and  outputs  as  a  column  vector  by  stacking 
successive  columns  of  the  matrix  below  the  previous  column  as  shown  in  Eq.  (47)  (Hace 
&  Johnson,  1991): 
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(47) 


a' 

b 
c 

cL. 

Complete  Angle  Scatter  Instrument  (CASI) 

AFIT's  complete  angle  scatter  instrument  (CASI)  is  used  to  experimentally 
determine  the  BRDF  of  selected  materials.  A  calibrated  laser  illuminates  the  sample  at 
different  incident  angles,  and  the  resultant  reflections  from  the  sample  {and/or  the 
transmission  through  the  sample)  are  measured  and  recorded  by  a  sensor  mounted  on  a 
movable  arm.  Both  in-plane  and  out-of-plane  measurements  can  be  performed  based  on 
the  geometry  of  the  sample's  orientation.  Figure  10  shows  AFIT's  CASI. 


Figure  10.  AFIT's  complete  angle  scatter  instrument. 
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Conclusion 


This  chapter  provided  a  description  of  dual  photography  as  well  as  a  review  of  the 
radiometry,  the  BRDF,  the  MTF  and  linear  algebra.  The  next  chapter  uses  the  concepts 
discussed  under  the  radiometry  section  to  develop  the  mathematical  theory  of  dual  and 
indirect  photography.  Chapter  IV  uses  the  concepts  of  the  modified  MTF  to  evaluate  the 
image  quality  of  the  images  produced  using  the  dual  or  indirect  techniques.  Chapter  V 
uses  the  linear  algebra  to  detail  a  possible  method  to  increase  image  quality  over 
traditional  blind  deconvolution  techniques. 
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III.  Radiometric  Theory  of  Indirect  Photography 


This  chapter  presents  the  radiometric  theory  of  indirect  photography.  First,  the 
radiometric  theory  of  dual  photography  is  developed;  the  theory  is  then  expanded  to  the 
radiometric  theory  of  indirect  photography.  The  chapter  closes  out  with  a  small  angle 
approximation  of  dual  photography.  The  contents  of  this  chapter  dealing  with  the  theory 
of  dual  and  indirect  photography  were  originally  presented  at  the  SPIE  Conference  on 
Reflection,  Scattering,  and  Diffraction  from  Surfaces  II  (Hoelscher  &  Marciniak,  2010). 
It  has  also  been  submitted  for  publication  to  Optics  Express  (Hoelscher  &  Marciniak, 
2011). 

Dual  Photography 

To  radiometrically  model  the  dual  photography  experiment  shown  in  Figure  1,  the 
setup  in  Figure  1 1  will  be  used.  In  this  configuration,  a  laser  is  used  as  the  illumination 
source  instead  of  a  pixilated  projector.  Additionally,  instead  holding  the  object  fixed  and 
moving  the  laser  spot  and  camera  in  unison,  the  laser  and  camera  are  fixed  and  the 
playing  card  is  translated.  Four  coordinate  systems,  three  fixed  with  respect  to  each  other 
and  one  fixed  to  the  object,  are  used.  The  x  coordinate  system  is  a  fixed  coordinate 
system  in-plane  with  the  object's  translation.  The  laser  spot  is  centered  at  its  origin.  This 
coordinate  system  will  be  referred  to  as  the  fixed  object  frame  of  reference.  The  x' 
coordinate  system  is  attached  to  the  object,  i.e.  the  playing  card  in  the  original 
experiment,  with  the  center  of  the  object  at  the  origin.  This  is  the  only  coordinate  system 
that  changes  with  respect  to  any  other  coordinate  system  during  the  creation  of  the  dual 
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image  and  is  referred  to  as  the  object  frame  of  reference.  The  y  coordinate  system  is 
attached  to  the  non-specular  surface  imaged  by  the  camera,  i.e.  the  book  in  the  original 
experiment,  and  will  be  referred  to  as  the  imaged  reflector.  The  z  coordinate  system  is 
fixed  and  attached  to  the  lens  of  the  imaging  system. 


Figure  11.  Dual  photography  coordinate  systems 


Using  the  setup  in  Figure  11,  the  irradiance  due  to  the  laser  in  the  fixed  object 
frame  of  reference  in  its  most  general  form  is  Eobj{dobj,(pobj,x,?i)  where  dobj  is 
incident  elevation  angle  of  the  irradiance  with  respect  to  the  surface  normal,  <pobj  is  the 
incident  azimuth  angle  of  the  irradiance  with  respect  to  the  fixed  object  frame  of 
reference,  x  is  a  2-D  vector  which  denotes  the  position  in  the  fixed  object  frame  of 
reference  and  A  is  the  wavelength  of  the  irradiance.  Given  the  irradiance  in  the  object 
frame  of  reference  and  using  Eq.  (26),  the  radiance  from  the  object  frame  of  reference 
can  be  written  as: 
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(48) 


Lobj  (@obj>  tpobj*  @  objf  0  Q^jt  ~ 

^obj^objf  tpobj’  tyfobj  (p  obj>  tpobj >  9  obj>  4*  obj>  X,  X  ,  A 

where  9'0bj  is  the  reflected  elevation  angle  with  respect  to  the  normal  of  the  fixed  object 
frame  of  reference,  4>'obj  is  the  reflected  azimuth  angle  with  respect  to  in  the  fixed  object 

frame  of  reference,  x'  is  the  offset  between  the  fixed  object  frame  of  reference  and  the 
object,  and  fobJ-  is  the  BRDF  of  the  object.  Because  the  irradiation  source  is  a  laser,  the 
wavelength,  A,  will  be  considered  constant  and  dropped  from  further  equations  for 
brevity.  Additionally,  because  a  laser  is  used  as  the  irradiation  source,  the  incident  solid 
angle  can  be  considered  constant  and  the  irradiance  in  the  fixed  object  frame  of  reference, 
Eobj,  can  be  written  solely  as  a  function  of  position  in  the  fixed  object  frame  of  reference, 
x.  Therefore,  Eq.  (48)  can  be  simplified  to: 


L()b  j  ("T  X  ,  9obj,  (pobj’  ^  obj)  ~  Eobj  Oc)fobj(x>  %  >  @objf  tpobjt  n'obj.)  (49) 

where  the  reflected  elevation  and  azimuth  angles,  9'0bj  and  <p'obj,  have  been  combined  to 
form  the  reflected  solid  angle,  Cl'obj.  Assuming  the  BRDF  of  the  object  is  isotropic  and 
uniformly  scaled  in  magnitude  by  the  reflectance  of  the  object  at  that  point  two 
simplifications  can  be  made:  (1)  the  BRDF  of  the  object  no  longer  has  a  dependence  on 
the  incident  azimuth  angle  with,  <pobj  and  (2)  the  BRDF  can  be  decomposed  into: 


fobjijX’X’^obj’^obj')  P(.X  X^fpb[x,Qobj,£lobj') 


where  fph  is  the  underlying  angular  shape  of  the  BRDF  (sometimes  referred  to  as  the 
"phase  function"  and  therefore  the  ph  subscript)  that  is  scaled  by  p,  the  reflectance  of  the 
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object  at  the  point,  x' .  To  validate  the  assumption  that  the  BRDF  can  be  decomposed 
into  a  reflectance  and  a  phase  function,  AFIT’s  CASI  was  used  to  measure  the  BRDF  of 
the  white,  black  and  red  portions  of  a  standard  playing  card,  the  results  can  be  found  in 
Appendix  A.  Substituting  Eq.  (50)  into  Eq.  (49)  the  radiance  of  the  object  becomes: 

Lobj(%,X  i^obj’^obj')  EobjCx^p^X  %')fph(^obj’^‘obj') 


Given  the  radiance  of  the  object  reflector,  and  using  Eq.  (19),  the  differential  irradiance 
on  the  imaged  reflector  from  a  differential  area  on  the  object  is: 


dEim  (x,  X  ,  0ojjj, 

^objf^irn)  Lobji^X,  X  ,Q0bj >  ^ ob j )  COS  dim  dElim 


(52) 


where  0im  is  the  incident  elevation  angle  with  respect  to  the  surface  normal  of  the 
irradiance  on  the  imaged  reflector  and  dQ.im  is  the  differential  solid  angle  incident  on  the 
imaged  reflector  that  is  subtended  by  a  differential  projected  area  of  the  object  (see 
Figure  12). 


Figure  12.  Differential  solid  angle  ( dSlim )  and  incident  angle  with  respect  to  the 

normal  ( 0im ) 
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Given  the  fixed  geometry  of  the  experiment,  Q.obj,  9im  and  dQ.im  depend  only  on 
the  positions  in  the  fixed  object  frame  of  reference  and  on  the  imaged  reflector,  therefore 
they  can  be  written  solely  as  functions  of  x  and  y  (see  Figure  13). 


Figure  13.  Angular  dependence  on  x  and  y. 


Rewriting  Eq.  (52)  explicitly  in  terms  of  x  and  y: 

dEim{x,x,y,eobj)  =  Lobj(x,x,y,6obj)cos6im(x,y)  dQim(x,y ) 


(53) 


Using  the  definition  of  a  solid  angle,  Eq.  (7),  the  differential  solid  angle  incident  on  the 
imaged  reflector,  dQ.im,  can  be  rewritten: 


cos  0ohi(x,  y) 

dQim(x,y)  = - 2  {  dhobj(x) 

r4(**y) 


(54) 


where  90bj  is  the  reflected  elevation  angle  of  the  radiance  with  respect  to  the  surface 

normal  of  the  differential  area  of  the  object,  rim  is  the  range  between  points  in  the  fixed 
object  frame  of  reference  and  the  imaged  reflector,  x  and  y,  respectively,  and  dAobj  is 
the  differential  area  of  the  object  (see  Figure  14). 
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X/ 


Figure  14.  Converting  from  a  differential  solid  angle  to  a  differential  area 


Substituting  Eq.  (54)  into  Eq.  (53)  yields: 
dEim(x,x  ,y,  d0bj) 

r-  -  n  \  n  { cose'obj(x,y) 

L0bj\.x>x  ,y, 6obj)cosOim(x, y)  [  dAobj{x) 


(55) 


By  combining  like  terms,  Eq.  (55)  can  be  rewritten: 

dEim{x,x,y,90b] )  =  Lobj{x,x,y,0obj)a(x,y)dAobj(x ) 


(56) 


where 


a  Or,  y) 


cos0ob;-  (x,  y)cosOim  (x,  y) 


(57) 


Again,  given  the  irradiance  on  the  imaged  reflector  and  using  Eq.  (26),  the  radiance  from 
the  imaged  reflector  can  be  written  as: 


,y.Z,  Qobj)  Eim{x,x  ,  y,  dobj)  fim(y>  ^irrv  <Plm>  dim’  <Pim) 


(58) 


where  fim  is  the  BRDF  of  the  imaged  reflector,  9im  and  Q'im  are  the  incident  and  reflected 
elevation  angles  with  respect  to  the  normal  of  the  imaged  reflector,  <pim  and  (p'im  are  the 
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incident  and  reflected  azimuth  angles  with  respect  to  the  imaged  reflector  frame  of 


reference,  and  y  is  2-D  vector  that  denotes  the  position  on  the  imaged  reflector. 
Assuming  the  imaged  reflector  is  uniform  and  isotropic,  the  BRDF  can  be  simplified  and 
Eq.(58)  can  be  rewritten: 


LimipC’X  >y>z>  Qobj)  Eim(%>X  >y >@obj')  fim(.@im>  ^-im) 


(59) 


where  the  reflected  elevation  and  azimuth  angles,  Q[m  and  (p[m ,  have  been  combined  to 
form  the  reflected  solid  angle  Q[m.  Converting  9im  and  fl[m  into  positions  in  the  object 
frame  of  reference,  x,  position  on  the  imaged  reflector,  y,  and  position  on  the  lens  of  the 
imaging  system,  z,  Eq.  (59)  can  be  rewritten: 

Lim  (x,  x y,  z,  Gobj )  =  Eim  (x,  x y,  Gobj )  fim  (x,  y,  z )  (60) 


Given  the  radiance  from  the  imaged  reflector,  and  again  using  Eq.  (19),  the  differential 
irradiance  at  any  point  on  the  lens  of  the  imaging  system  from  a  differential  area  on  the 
imaged  reflector  is: 


dElens(x,x , y, z,  9obj  )  =  Lim{x,x' ,y,z,  9 obj)cos9lens(y,z)dQ.lens(y,z) 


(61) 


where  9iens  is  the  incident  elevation  angle  with  respect  to  the  surface  normal  of  the  lens 
of  the  imaging  system  and  d£Llens  is  the  incident  solid  angle  of  the  radiation  on  the  lens 
subtended  by  a  differential  area  of  the  imaged  reflector.  Using  Eq.  (7)  to  convert  the 
solid  angle,  Eq.  (61)  can  be  rewritten: 
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dEiens(^X,  X  >y >  Z,  @obj  ) 


, _ , _  .  _ /  cosdim(y,z )  '  ^  ^ 

Urn  (*,  X,y,z,  Gob  j )  cos9lens  (y,  z)  (  2  dAlm  (y) 

rlensyy>  ZJ 


where  9im  is  the  reflected  elevation  angle  with  respect  to  the  surface  normal  of  the 
imaged  reflector,  rlens  is  the  distance  between  the  points  on  the  imaged  reflector  and  the 
lens  of  the  imaging  system,  y  and  z  respectively,  and  dAim  the  is  differential  area  on  the 
imaged  reflector.  Again  combining  like  terms,  Eq.  (62)  can  be  simplified: 


d-Eiens  (x , x , y , z ,  0 obj  ^  Lim{x,x  ,y,  z,  90bj)P(y> z)dAim(y) 


(63) 


where 


Piy.z) 


cos0-m(y,z)cos6>tens(y,z) 

r?ens(y>z) 


(64) 


Given  Eq.  (63)  and  by  using  Eqs.  (51),  (56),  and  (60)  the  irradiance  on  a  point  on 
the  lens  of  the  imaging  system,  z,  for  a  given  object  position  and  incident  elevation  angle 
of  the  laser  on  the  fixed  object  frame  of  reference,  x  and  6obj,  respectively,  is: 


II 

im  obj 


EobjMp (x'  -  x)fph{x,  y,  d0bj)fim(x,  y,  z)a(x,  y)  / 3{y,z)dx  dy 


(65) 


where  x  is  integrated  over  the  object,  and  y  is  integrated  over  the  imaged  reflector.  If  the 
irradiance  outside  the  laser  spot  is  zero,  Eq.  (65)  can  be  rewritten: 
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(66) 


E0bj(.x)p(x  -  x)fph(x>y> dobj)fim(x, y, z)a(x, y)  p(y,z)dx  dy 

im  las 

where  x  is  now  integrated  over  the  laser  spot.  Given  Eq.  (66),  the  total  flux  collected  by 
the  lens  for  a  given  pixel  on  the  camera  can  be  written  as: 


^pixj  (x  >  @ob j)  ~ 


I  j  J  EobJ(x)p(x-x)fph(x,y,0obj)fim(x,y,zM^ 


(67) 


lens  fovi  las 

where  fovi  is  the  projected  area  of  camera  pixel  i  on  the  imaged  reflector  (see  Figure 
15). 


Figure  15.  Field  of  view  of  an  individual  pixel 


The  order  of  integration  can  be  rearranged  and  Eq.  (67)  simplified  to: 


where: 


rj (x,0obj)=  I  Eobj(x)fph(x,y,eobj)fim(x,y,z)a(x,y)  p(y,z)dy  dz  (69) 
Hens  J fovi 

As  a  consequence  of  Eq.  (68),  a  dual  image  can  be  created  by  using  any  single  pixel, 
group  of  pixels  or  the  entire  digital  image  without  explicit  knowledge  of  the  geometry,  as 
long  as  the  same  pixel  or  set  of  pixels  is  used  to  create  the  dual  image  across  all  of  the 
recorded  data  images.  Furthermore,  with  a  change  of  variables,  (x  =  x  —  x  ),  Eq.  (68) 
can  be  rewritten: 

®piXi(x  >  @obj)  = 

showing  H  to  be  the  convolution  kernel,  i.e.  the  point  spread  function  (PSF),  for  the  dual 
image.  If  the  irradiance  of  the  laser  spot  and  some  knowledge  of  the  BRDFs  and 
geometries  in  Tj  are  known,  the  quality  of  the  dual  image  can  be  improved  by  a 
deconvolution  of  the  dual  image  and  this  kernel.  It  is  this  improvement  in  the  image 
quality  by  the  deconvolution  of  the  irradiance  on  the  object  of  interest  that  will  be  shown 
can  be  exploited  to  expand  the  concept  of  dual  photography  into  one  of  indirect 
photography  and  allow  for  the  recovery  of  information  that  is  not  directly  visible  to  either 
the  controlling  illumination  source  or  the  imaging  system. 

Indirect  Photography 

As  previously  stated,  one  limiting  factor  of  dual  photography  is  the  requirement 
for  the  illumination  source,  or  the  imaging  system,  to  have  a  direct  view  of  the  object  of 


I  Ij(x’  -  x" ,  eobj )  p(x")  dx" 
Has 


(70) 
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interest.  To  eliminate  that  requirement,  the  dual  photography  experiment  modeled  above 
will  now  be  expanded  and  modeled  with  the  illumination  source  reflecting  from  a  non- 
specular  surface  prior  to  illuminating  the  object  of  interest.  For  example,  if  the  laser  is 
moved  adjacent  to  the  camera  so  that  it  could  not  illuminate  the  face  of  the  playing  card 
directly  but  now  illuminates  the  diffuse  imaged  reflector  as  shown  in  Figure  16. 


projector 


Figure  16.  Real  world  setup  for  (a)  dual  photography  and  (b)  indirect  photography. 


To  aid  in  the  modeling,  the  first  surface  has  been  separated  from  the  imaged 
reflector  and  an  additional  fixed  reference  frame,  w,  has  been  added  to  describe  the  first 
non-specular  surface  and  will  be  referred  to  as  the  wall  reflector  (see  Figure  17). 

As  in  the  previous  section,  because  the  illumination  source  is  a  laser,  both  the 
wavelength  and  incident  solid  angle  can  be  considered  constant  and  the  general  form  of 
the  irradiance  on  the  wall  reflector  EwaU  (9waU,  <pwau,  w,  X)  can  again  be  simplified  to 
Ewau(w),  where  w  is  a  2-D  vector  which  denotes  the  position  on  the  wall  reflector. 
Given  the  irradiance  on  the  reflector  and  using  Eq.  (26),  the  radiance  from  the  wall 
reflector  is: 
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Lwalliyv ,  9wall>  (pwalb  9  wall’  0  wall )  — 

(71) 

Ewall(w)fwaii{w,  @wall’  tpwall'  9  wall’  0  wall ) 

where  0wa;;  and  9'wall  are  the  incident  and  reflected  elevation  angles  with  respect  to  the 
surface  normal  of  the  wall  reflector,  (pwaii  and  <p'waii  are  the  incident  and  reflected 
azimuth  angles  with  respect  to  the  wall  reflector  frame  of  reference  and  fwaU  is  the 
BRDF  of  the  wall  reflector. 


Figure  17.  Indirect  Photography  coordinate  systems 


Assuming  the  wall  reflector  is  both  homogenous  and  isotropic,  the  BRDF  of  the 
wall  reflector  can  be  simplified  and  Eq.  (71)  can  be  rewritten: 

Lwall(yv >  ^wall’  ^  wall)  ^wa!I (^)/wa!i (^wa!i>  ^  wall)  (72) 


where  the  reflected  elevation  and  azimuth  angles,  9'wan  and  0'  have  been  combined 
to  form  the  reflected  solid  angle,  £l'waii-  Given  the  radiance  from  a  differential  area  of 
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the  wall  reflector  and  using  Eq.  (19),  the  differential  irradiance  on  the  object  frame  of 
reference  is: 

dEobj(w,  Qwall>  Oobj’  &obj)  Ewau(w,  9Wall>  0.  wall)COS9objdQ.obj  (73) 


where  90bj  is  the  incident  elevation  angle  of  the  irradiance  with  respect  to  the  surface 
normal  in  the  fixed  object  frame  of  reference  and  dCLobj  is  the  differential  incident  solid 
angle  of  the  irradiance.  Again,  using  Eq.  (7)  and  converting  the  differential  solid  angle  to 
differential  area  yields: 


dEob  j  (w,  dob  j,  9waU,  Elwall)  — 


Lwall(yv ,  9 wall,  fl  Wall)COSd obj 


^COS  9, 


(74) 


wall 


dA 


wall 


'obj 


where  9waU  is  the  reflected  elevation  angle  with  respect  to  the  surface  normal  of  the 
radiance  from  the  wall  reflector,  robj  is  the  range  between  the  points  on  the  wall  reflector 
and  the  position  in  the  fixed  object  frame  of  reference,  w  and  x  respectively,  and  dkwaa 
is  the  differential  area  of  the  wall  reflector.  As  in  the  previous  section,  because  of  the 
fixed  geometry  between  the  wall  reflector  and  the  fixed  object  frame  of  reference,  all  of 
the  angles  can  be  written  explicitly  as  functions  of  w  and  x.  Rewriting  Eq.  (74)  yields: 


dEob  j (w ,  X,  9 wall) 

Ewall  (w,  X,  9 wall )  COs9obj  (w,  x) 


COs6>wct;;(w,x) 

robj(™’  *) 


(75) 


dAWaii(yv) 
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and  by  combining  like  terms,  Eq.  (75)  can  be  rewritten: 


dE0bj(w,  x,  6waii)  —  Lwau  (w,  x,  ewaU)  y(w,  x )  dAwall(w) 


(76) 


where: 


y(w,  x ) 


cos  6>wct„  (w,  x)cos0ob;-  (w,  x) 
r^jiw.x) 


(77) 


Eq.  (76)  can  now  be  inserted  into  Eq.  (67)  to  give  the  total  flux  collected  by  a  single  pixel 
for  a  given  object  position  and  the  incident  elevation  angle  of  the  laser  with  respect  to  the 
normal  of  the  wall  reflector: 


’  @wall)  ~  I  I  I  I  Ewap(w)  p(x  X^  fwau(w,  X,  9wan) 
Has  J obj  Jfovi  Hens 


(78) 


fph(w,x,y)fim(x,y,z)y(w,x)a(x,y)  0(y,z)  dz  dy  dx  dw 


where  dz  is  integrated  over  the  area  of  the  lens  of  the  imaging  system,  dy  is  integrated 
over  the  projected  area  of  camera  pixel  i  on  the  imaged  reflector,  dx  is  now  integrated 
over  the  entire  object  and  dw  is  integrated  over  the  laser  spot.  Eq.  (78)  can  be  simplified 
to: 


d> 


pixt 


(.X1,  ewall )  =  I  I  Ii(w,x)  p(x'  -  x)  X 
Has  Jobj 


(79) 


Ewaii(w)fwaii(W’X,  dwall)y(w,x)  dx  dw 


where: 


fi(w,x) 


-I  s 

j  fnv:  J  h 


fovi  Jlens 


fph  O,  x,  y)fim  (x,  y,  z)  a(x,  y)  /3(y,z)  dz  dy 


(80) 
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The  order  of  integration  can  be  rearranged  and  Eq.  (79)  further  simplified  to: 


^pixjx  ’  @ wall ) 


=  /  T> 

Jobj 


(.X,  6waii)  p(x'  -x)  dx 


(81) 


where: 


T;(x,  ewaU)  =  I  T^w.x,  ewall)  Ewall(w)fwall(w,x,  ewall )  y(w,x)dw 

Jlas 


(82) 


As  with  Eq.  (68),  a  consequence  of  Eq.  (81)  is  an  indirect  image  can  be  created  by 
using  any  single  pixel,  group  of  pixels  or  the  entire  digital  image  without  explicit 
knowledge  of  the  geometry,  as  long  as  the  same  set  of  pixels  is  used  to  create  the  dual 
image  across  all  of  the  recorded  images.  Likewise,  Eq.  (82)  can  also  be  rewritten  with  a 
change  of  variables: 

^pix i(.X  >  dwall )  =  f  Tj(x  —  X  ,  dwan )  p(x  )  dx  (^3) 

Jobj 


as  in  the  case  of  the  dual  image,  the  indirect  image  is  the  convolution  of  the  point  spread 
function,  Ti(  and  the  object  of  interest,  p.  Due  to  the  nature  of  the  problem,  the  point 
spread  function  will  probably  not  be  fully  known;  it  should  however,  be  possible  to 
improve  the  image  quality  using  blind  deconvolution  techniques. 


Dual  Photography  Approximation 

In  the  above  derivations  of  dual  and  indirect  photography,  assumptions  were 
made  about  the  BRDFs  of  the  reflectors,  i.e  that  they  were  homogeneous  and  isotropic. 
No  assumptions,  however,  were  made  about  the  geometry  of  the  setup(s).  If  two 
assumptions  are  made  about  the  geometry  of  the  setup,  Eq.  (70)  can  be  significantly 
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simplified  and  the  convolution  kernel  in  dual  photography  can  be  approximated  with  the 
laser  irradiance  on  the  object. 

The  two  assumptions  that  must  be  made  are:  (1)  both  the  area  of  the  laser  spot  on 
the  object  and  the  projected  area  of  pixel  i  on  the  imaged  reflector  are  much  less  than  the 
range  between  the  two  points  squared,  rfm,  and  (2)  the  incident  elevation  angle  of  the 
laser  on  the  object,  0obj,  is  constant.  The  first  assumption  means  rim  can  be  considered 
constant  for  the  range  between  any  point  with  the  laser  spot  and  any  point  in  within  the 
projected  pixel  area.  Also  if  the  first  assumption  is  true,  then  the  reflected  elevation 
angle  from  object  frame  of  reference,  O'obj,  will  vary  only  slightly  across  all  possible 
combinations  of  x  and  y;  therefore,  cos6'obj(x,y),  of  the  angle  can  be  considered 
constant.  Likewise,  the  cosine  of  the  incident  elevation  angle  on  the  imaged  reflector, 
9im,  can  also  be  considered  constant.  This  allows  Eqs.  (57)  and  (64)  to  be  rewritten: 


a(x,y) 


cosffipfr  j  (Xp,  y0)cosfl,m(x0,  y0) 


(84) 


P(y,z) 


cos6>-m(y0,  z)cQ50;en5(yo,  z) 

rl2ens(y0.z) 


=>&(*) 


(85) 


where  x0  and  y0  are  constants  representing  the  center  position  of  the  laser  spot  and  the 
pixel,  respectively,  and  the  i  subscript  denotes  the  pixel.  If  the  second  assumption  is  true, 
as  is  the  case  when  the  laser  is  translated  horizontally  and  vertically,  and  by  holding  x 
and  y  constant,  the  BRDF  phase  component  of  the  object,  fph,  can  also  be  considered 
constant;  i.e.  if  the  reflected  solid  angle  is  small,  the  BRDF  within  that  solid  angle  will 
vary  only  slightly,  and  can  therefore  be  rewritten: 
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(86) 


fph(x>  y>  Oobj)  ^  fph  y^Ofyof  @0 obj)  ^  fp^i 

Again,  if  the  laser  spot  is  small  in  comparison  to  the  range  between  the  laser  spot  and  the 
projection  of  the  pixel,  then  the  incident  solid  angle  of  the  irradiance  on  the  imaged 
reflector  can  be  considered  constant  across  the  laser  spot  and  the  BRDF  of  the  imaged 
reflector,  fim,  can  be  rewritten: 

fim(x,y,z)  =>  fim(x0,y0,z )  =>  fim.{z )  (87) 

Using  Eqs.  (84)  through  (87),  Eq.  (69)  becomes: 


F'j(x)  =  Eobj(x)  fph.fim.(z)  ai  pi(z)dy  dz  (88) 

Jlens  J  fovi 

Rearranging  the  integral  yields: 

E'i(x)  =  Eobj(x)aifph.  f  dy  (  firn.(z)Pi(z)dz  (89) 

Jfovi  Jlens 


The  integration  over  the  field  of  view  of  the  pixel  can  be  evaluated  and  yields  the  area  of 
the  pixel  which  is  a  constant. 


and  given  the  fixed  geometry  of  the  setup,  the  integration  of  fim.  and  over  the  area  of 
the  lens  can  be  evaluated  and  yields  a  constant 
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Substituting  Eqs.  (90)  and  (91)  into  Eq.  (89)  yields: 

r'i(x)  =  Eobj(x)a.i  fphAifirn.  (92) 

which  can  be  simplified  to: 

r'i(x)  =  Kt  Eobj(x)  (93) 

where 

K,  =  a,  fphiA,f,mi  <94) 

It  is  important  to  note  that  Kt  varies  from  one  pixel  to  the  next  but  for  each  pixel  i,  Kt  will 
remain  constant  for  all  data  images.  Substituting  Eq.  (94)  into  Eq.  (68)  yields: 

)  = 

which  through  a  change  of  variables  can  be  rewritten: 

=  Ki  [  Eobj{x'  -  x)p{x")dx  (96) 

Has 

Eq.  (96)  suggests  the  complete  geometry  of  a  dual  photography  setup  is  not  necessary  to 
improve  the  image.  Instead,  the  irradiance  of  the  controlling  illumination,  in  this  case  the 
laser,  can  be  used  as  the  deconvolution  kernel. 

Conclusion 

This  chapter  developed  the  radiometric  theory  of  both  dual  and  indirect 
photography.  It  also  included  a  simplification  of  the  radiometric  equation  of  dual 
photography.  To  verify  Eqs.  (70),  (83)  and  (96),  i.e.  that  the  image  produced  by  either 
dual  or  indirect  photography  is  a  convolution  of  the  original  object  of  interest  and  either 


I  Eobj{x)p(x' —  x)dx 
Has 


(95) 
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the  laser  spot  in  the  case  of  dual  photography  or  an  unknown  point  spread  function  in  the 
case  of  indirect  photography,  both  a  MATLAB  simulation  and  a  physical  experiment 
described  in  the  next  chapter  of  this  document  were  used. 
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IV.  Experimental  Verification 


To  verify  the  radiometric  theory  of  indirect  photography  described  in  Chapter  III, 
a  basic  simulation  and  three  different  experiments  were  accomplished  in  a  building  block 
approach.  The  first  experiment  was  a  1  -D  experiment  with  the  illumination  source  and 
the  camera  separated  to  mimic  the  theoretical  model.  Sinusoidal  slides  and  a  separated 
object  reflector  were  used  as  the  object.  The  experiment  was  then  expanded  to  2-D 
where  the  sinusoidal  slides  were  replaced  with  more  representative  2-D  objects.  In  the 
final  experiment,  the  illumination  source  and  the  laser  were  co-located  to  create  a  real- 
world  scenario.  This  chapter  describes  the  simulation  as  well  as  the  three  experimental 
setups  and  discusses  the  results. 

Simulation 

Following  the  development  of  the  radiometric  theory  of  indirect  photography,  a 
MATLAB  simulation  was  created  to  verify  the  results  of  Eqs.  (70),  (83)  and  (96),  i.e.  that 
the  image  quality  of  both  dual  and  indirect  images  are  improved  following  the 
deconvolution  process.  To  that  end,  Eqs.  (67)  and  (78)  were  used  to  calculate  the  total 
flux  incident  on  the  lens  of  the  imaging  system  for  the  dual  and  indirect  image, 
respectively.  The  resultant  images  were  then  improved  through  the  deconvolution 
processes.  The  image  quality  was  evaluated  for  both  the  unimproved  images  and 
recovered  images  following  deconvolution.  A  description  of  the  simulation  and  the 
subsequent  results  are  described  below. 
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The  simulation  was  based  on  the  micro-facet  model  discussed  in  the  background 
section  of  this  document  and  consisted  of  a  wire  frame  model  for  each  reflector.  The 
three  reflectors  and  the  lens  of  the  imaging  system  were  modeled  by  the  following  sets  of 
points:  (1)  the  wall  reflector  consisted  of  9  points  in  a  3x3  pattern,  (2)  the  object  reflector 
consisted  of  625  points  in  a  25x25  pattern,  (3)  the  imaged  reflector  consisted  of  441 
points  in  a  21x21  pattern,  and  (4)  the  lens  of  imaging  system  consisted  of  9  points  in  a 
3x3  pattern.  Two  objects  were  used  for  the  simulation.  The  first  object  was  a  single 
white  square  on  a  black  background,  so  the  PSF  of  the  system  could  be  evaluated.  The 
second  object  consisted  of  a  two  white  squares  separated  both  vertically  and  horizontally 
by  one  side  length  on  a  black  background.  For  reference,  the  squares  used  to  create  the 
objects  subtended  approximately  0.08  mrad  at  the  distance  modeled  in  the  simulation. 

Each  component  of  Eqs.  (67)  and  (78)  was  calculated  and  stored  in  look-up  tables 
to  decrease  the  time  requirements  to  run  the  model.  The  Priest  and  Meier  BRDF, 
Eq.  (27),  a  well  studied  form  of  the  BRDF  using  the  glint  vector,  was  used  to  used  to 
model  the  BRDF  of  each  of  the  reflectors.  The  BRDF  of  the  wall  reflector  was  chosen  so 
that  when  the  object  was  in  the  center  position,  the  radiance  from  the  wall  reflection 
would  cover  the  entire  object.  Figure  18  shows  the  glint  angle  from  the  center  point  on 
the  wall  reflector  to  every  point  in  the  fixed  object  frame  of  reference  and  Figure  19  is  an 
overlay  of  the  irradiance  in  the  object  frame  of  reference  and  the  two  square  object  when 
it  is  in  the  center  position. 

The  laser  was  modeled  as  a  3x3  Gaussian  beam,  Eq.  (97),  shows  the  matrix 
representation  of  the  irradiance  on  the  wall  reflector  used  for  the  simulation. 
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Figure  18.  Glint  Angle  from  the  center  of  the  wall  reflector  to  every  point  on  the 

object  reflector 


Figure  19.  Simulated  object  irradiance 
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The  object  was  then  placed  in  the  upper  left  comer  of  the  fixed  object  frame  of 
reference  and  beginning  from  the  wall  reflector,  every  possible  path  to  the  imaging 
system  was  evaluated,  i.e.  from  each  point  on  the  wall  reflector  to  every  point  on  the 
object  reflector  to  every  point  on  the  imaged  reflector,  etc.  The  total  flux  impinging  on 
the  lens  of  the  imaging  system  was  summed  to  simulate  the  entire  image  was  being  used 
to  create  the  dual  or  indirect  image.  The  object  was  then  translated  horizontally  and 
vertically  through  each  of  the  possible  positions  and  the  process  repeated.  The  dual  and 
indirect  images  created  by  the  simulation  are  shown  in  Figure  20.  The  full  MATLAB 
code  used  to  create  the  images  can  be  found  in  Appendix  B. 


(a) 


\ 


(c)  (d) 

Figure  20.  Simulation  results  of  the  single  square  (a)  dual  and  (b)  indirect  images 
and  the  two  square  (c)  dual  and  (d)  indirect  images 


As  stated  earlier,  Eqs.  (70),  (83)  and  (96)  suggest  that  both  dual  and  indirect 
image  quality  can  be  improved  through  a  deconvolution  process.  The  dual  images  were 
improved  using  two  of  MATLAB's  deconvolution  algorithms,  deconvlucy  which  is  based 
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on  the  Lucy-Richardson  (LR)  (Richardson,  1972)  (Lucy,  1974)  method  and  a  blind 
deconvolution  algorithm  deconvblind  (Holmes,  1992)  (Gonzalez,  Woods,  &  Eddins, 
2004,  pp.  176,179).  The  irradiance  of  the  laser  shown  in  Eq.  (97)  was  used  as  the  PSF 
(deconvolution  kernel)  in  the  deconvlucy  routine.  Because  explicit  knowledge  of  the 
entire  setup  would  be  required  to  fully  develop  the  PSF  (deconvolution  kernel)  for 
indirect  photography,  which  is  unrealistic  in  a  real  world  scenario,  only  the  blind 
deconvolution  was  used  on  the  indirect  images.  Following  each  iteration  in  the 
deconvolution  process,  Eq.  (36)  was  used  to  quantify  image  quality  of  each  iteration's 
recovered  image.  The  algorithm  was  allowed  to  run  for  ten  iterations,  after  which  there 
was  negligible  improvement  with  each  successive  iteration.  Figure  21  shows  the  image 
quality  after  each  iteration  of  the  blind  deconvolution  algorithm  for  both  of  the  indirect 
images. 


Figure  21.  Simulation  image  quality  improvement  per  deconvolution  iteration 
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Table  4  summarizes  the  results  of  the  image  quality  assessment  in  the  original  and 
improved  images.  Figure  22  shows  the  recovered  images  after  the  blind  deconvolution 
process.  (Note:  Because  there  was  no  discemable  difference  between  the  recovered 
images  using  the  LR  or  blind  deconvolution  algorithm,  only  those  recovered  using  the 
blind  deconvolution  are  shown.) 


Table  4.  Simulation  image  quality2 


Raw  Image 

LR  Deconv 

Blind  Deconv 

1  point  dual  image 

0.8326 

0.9999 

0.9998 

1  point  indirect  image 

0.6616 

n/a 

0.8160 

2  point  dual  image 

0.8326 

0.9999 

0.9989 

2  point  indirect  image 

0.6863 

n/a 

0.8150 

(c)  (d) 

Figure  22.  Simulation  results  following  the  blind  deconvolution  of  the  single-square 
(a)  dual  and  (b)  indirect  images,  and  the  two-square  (c)  dual  and  (d)  indirect  images 


2  The  image  quality  calculations  were  carried  out  to  four  significant  digits  to  show  the  difference  between 
the  LR  and  the  blind  deconvolutions.  In  a  real  world  scenario,  the  noise  floor  will  be  the  determining  factor 
on  the  number  of  significant  digits  the  equations  will  support. 
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While  somewhat  limited  in  scope  due  to  computer  processing  speed  and  memory 
requirements,  the  above  results  verify  that  the  image  quality  of  both  dual  and  indirect 
images  as  modeled  by  Eqs.  (67)  and  (78)  can  be  improved  through  deconvolution.  It  also 
validated  the  small  angle-approximation  of  dual  photography  and  verified  Eq.  (96)  could 
be  used  to  improve  the  image  quality  of  dual  images  if  the  irradiance  in  the  object  frame 
of  reference  is  known.  Due  to  time  constraints,  a  full  scale/higher  fidelity  version  of  the 
simulation  was  not  pursued,  opting  instead  to  begin  a  physical  experiment.  The 
following  section  describes  the  1-D  experiment  and  its  results  in  detail. 

1-D  Experimental  Setup 

The  dual  and  indirect  photography  1-D  experiments  were  set  up  in  accordance 
with  Figure  23  (a)  and  (b),  respectively,  where  a  633nm  HeNe  laser  is  used  as  the 
illumination  source.  The  imaged  reflectors,  and  wall  reflector  for  the  indirect  images, 
were  polished  aluminum  plates  that  had  been  finished  in  one  of  three  ways:  (1)  spray- 
painted  semi-gloss  white  paint,  (2)  spray-painted  with  a  flat  white  paint  or  (3)  polished 
with  600-grit  sandpaper  and  left  unpainted.  The  BRDF  of  each  of  the  reflectors  was 
measured  using  AFIT’s  CASI  instrument  and  the  respective  measurements  can  be  found 
in  Appendix  C.  The  object  reflector,  which  provides  the  phase  function  dependence,  fph, 
was  also  a  polished  aluminum  plate  with  a  flat  white  finish  similar  to  finish  of  the  second 
set  of  reflectors. 

Six  sinusoidal  slides  with  spatial  frequencies  from  0.1  to  3.0  cycles/mm  (0.02  to 
0.6  cycles/mrad  with  respect  to  the  wall-reflector-to-slide-distance)  were  used  as  the 
object  to  create  the  dual  and  indirect  images.  (See  Table  5  for  a  complete  list  of  slides 
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used.)  A  computer-controlled  translation  stage  moved  the  slide  through  the  path  of  the 


illumination  to  create  the  (x’  —  x)  dependence  for  each  image.  The  complete  laboratory 
set  up  is  shown  in  Figure  24. 


Table  5. 1-D  sinusoidal  slides 


slide 

Spatial  Frequency 
cycles/mm  cycles/mrad 

Slide  1 

0.1 

0.02 

Slide  2 

0.2 

0.04 

Slide  3 

0.5 

0.10 

Slide  4 

1.0 

0.20 

Slide  5 

2.0 

0.40 

Slide  6 

3.0 

0.60 

1,000  data  images  were  taken  with  the  sinusoidal  slide  translated  horizontally 
0.1mm  between  each  image.  These  data  images  were  then  used  to  create  the  dual  or 
indirect  image,  depending  on  the  setup,  in  accordance  with  the  dual  photography 
algorithm  described  in  the  background  section  of  this  document.  Figure  25  shows  the  0. 1 
cycle/mm  (0.02  cycle/mrad)  slide  being  illuminated  by  the  reflection  from  the  semi-gloss 
wall  reflector,  i.e.  when  creating  an  indirect  image. 

Representative  images  recorded  by  the  camera  through  a  non-transmissive  and 
transmissive  portion  of  the  0.1  cycles/  mm  slide  (0.02  cycles/  mrad)  dual  image  are 
shown  in  Figure  26  (a)  and  (b),  respectively. 
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Imaged  Reflector  (y) 


Mirror  Laser 

(b) 


Figure  23.  1-D  (a)  dual  photography  and  (b)  indirect  photography  setup 
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Reflector 


Wall 

Reflector 


-Camera 


Laser 


Turning  Mirror  Slide  Object  Reflector 

Figure  24. 1-D  Laboratory  setup 


Figure  25.  Indirect  illumination  on  the  0.1  cycles/mm  (0.02  cycles/mrad) 

slide 
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(a)  (b) 


Figure  26.  Data  image  of  0.1  cycles/mm  slide  through  (a)  a  non-transmissive  and  (b) 

transmissive  portion  of  the  slide 

To  improve  the  signal -to-noise  ratio,  the  entire  image  was  used  to  create  each  of 
the  dual  or  indirect  images.  The  average  intensity  of  each  digital  image  was  used  as  the 
total  irradiance  on  the  lens  for  each  position,  x ,  as  represented  by  Eqs.  (67)  and  (77)  for 
the  dual  and  indirect  images,  respectively. 

Beam  View  Analyzer  by  COHERENT®  was  used  to  obtain  the  cross-sectional 
power  distribution  of  the  laser  beam  used  as  the  illumination  source.  This  analysis  was 
used  to  create  an  estimation  of  the  object  irradiance,  Eobj,  as  a  function  of  position  which 
was,  in-tum,  used  as  the  deconvolution  kernel  for  MATLAB's  LR  deconvolution 
algorithm.  The  dual  images  created  with  the  semi-gloss  reflector  were  then  improved 
using  the  LR  and  the  blind  deconvolution  algorithms.  The  results  of  these 
deconvolutions  for  the  0.1  cycles/mm  slide  are  shown  in  Figure  27  while  an  expanded 
view  of  the  1.0  cycles/mm  slide  dual  image  and  deconvolutions  are  shown  in  Figure  28. 
Based  on  the  high  quality  match  of  the  two  deconvolution  techniques,  and  for 
consistency,  further  analysis  will  be  done  with  only  the  blind  deconvolution  algorithm. 
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Following  the  above  analysis,  both  the  dual  and  indirect  images  were  improved 
using  MATLAB's  blind  deconvolution  algorithm.  After  the  creation  of  the  improved 
images,  Fourier  analysis  was  accomplished  on  all  of  the  images  (dual,  indirect  and 
improved)  to  assess  the  amount  of  energy  in  the  fundamental  frequency  of  each  of  the 
corresponding  slides.  Figure  29-31  show  the  following  images  for  the  semi-gloss,  flat 
white  and  unpainted  reflectors,  respectively,  all  for  the  0.1  cycles/mm  (0.02  cycles/mrad) 
(a)  dual  images;  (c)  indirect  images;  (e)  overlay  of  the  dual  and  deconvolved  images;  and 
(g)  overlay  of  the  indirect  and  deconvolved  images;  (b),  (d),  (f),  and  (h)  are  the  Fourier 
transforms  of  (a),  (c),  (e),  and  (g),  respectively.  The  complete  set  of  dual  and  indirect 
images,  as  well  as  the  corresponding  plots  of  the  Fourier  transforms,  can  be  found  in 
Appendices  D-F  for  the  semi-gloss,  flat  white  and  unpainted  reflectors. 
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Figure  27.  0.1  cycles/mm  slide  dual,  and  dual  with  LR,  and  blind  deconvolutions 
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Figure  28.  Expanded  1.0  cycles/mm  slide  dual,  and  dual  LR,  and  blind 

deconvolutions 

As  previously  stated  the  modulation  transfer  function  (MTF)  is  one  technique  of 
quantifying  an  imaging  system's  ability  to  transfer  frequency  content  from  the  object  of 
interest  to  the  final  image.  Therefore,  the  transfer  of  information,  i.e.  spatial  frequencies, 
for  the  dual,  indirect  and  deconvolved  images  can  be  quantified  by  the  MTF  described  by 
Eq.  (34).  For  the  purposes  of  this  analysis,  the  energy  corresponding  to  the  fundamental 
frequency  of  each  slide  was  normalized  to  the  lowest  frequency  dual  image  (0.1 
cycles/mm  or  0.02  cycles/mrad  slide)  and  then  used  to  create  the  MTF.  The  experimental 
MTF  for  the  semi-gloss,  flat  white  and  unpainted  reflectors  are  shown  in  Figure  32-34. 

The  1-D  experiment  was  set  up  with  three  goals  in  mind.  First  was  to  validate  the 
theoretical  model  of  indirect  photography  and  insure  spatial  information  could  be 
recovered  from  an  indirect  image  created  from  wall  and  imaged  reflectors  with  different 
reflection  characteristics.  Second  was  to  verify  Eqs.  (70)  and  (83)  could  be  used  to 
improve  the  image  quality  of  both  the  dual  and  indirect  images.  The  final  goal  was  verify 
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the  small-angle  approximation  of  dual  photography,  and  therefore,  Eq.  (96)  could  be  used 
to  improve  the  image  quality  of  the  dual  image. 

The  creation  of  the  dual/indirect  images  and  subsequent  deconvolutions  as  shown 
in  the  MTFs  of  Figure  32  -34  verify  that  the  frequency  content  of  the  slides,  both 
visually  and  quantitative,  could  be  recovered  and  that  the  deconvolution  would  improve 
the  image  quality.  The  small-angle  approximation  was  also  confirmed  by  a  direct 
comparison  of  a  blind  deconvolution  and  a  FR  deconvolution  using  the  laser  profile  as 
the  kernel  for  the  deconvolution.  Having  verified  the  general  assumptions  and  validated 
the  radiometric  theory  in  1-D,  the  next  set  of  experiments  expanded  the  experiment  to 
2-D. 

2-Dimensional  Setup 

With  the  final  goal  of  a  real  world  setup  with  the  illumination  source  and  camera 
co-located  to  image  a  2-D  object,  an  intermediary  step  of  a  2-D  object  similar  to  the  1-D 
experiment  was  conducted.  The  sinusoidal  slides  and  object  reflector  were  replaced  by  a 
2-D  object  in  the  place  of  the  object  reflector  (see  Figure  35). 
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(g)  (h) 


Figure  29.  Semi-gloss  reflector  0.1  cycles/mm  (a)  dual  image  and  (b)  Fourier 
transform  (c)  indirect  image,  (d)  Fourier  transform  (e)  overlay  of  1.0  cycles/mm 
dual  and  improved  image  and  (f)  Fourier  transform  (g)  overlay  0.1  cycles/mm 
indirect  and  improved  image  and  (h)  Fourier  transforms 
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Figure  30.  Flat  white  reflector  0.1  cycles/mm  (a)  dual  image  and  (b)  Fourier 
transform  (c)  indirect  image,  (d)  Fourier  transform  (e)  overlay  of  1.0  cycles/mm 
dual  and  improved  image  and  (f)  Fourier  transform  (g)  overlay  0.1  cycles/mm 
indirect  and  improved  image  and  (h)  Fourier  transforms 
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Figure  31.  Unpainted  reflector  0.1  cycles/mm  (a)  dual  image  and  (b)  Fourier 
transform  (c)  indirect  image,  (d)  Fourier  transform  (e)  overlay  of  1.0  cycles/mm 
dual  and  improved  image  and  (f)  Fourier  transform  (g)  overlay  0.1  cycles/mm 
indirect  and  improved  image  and  (h)  Fourier  transforms 
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Figure  32.  Semi-gloss  reflector  1-D  experimental  MTF 
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Figure  33.  Flat  white  reflector  1-D  experimental  MTF 
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Figure  34.  Unpainted  reflector  1-D  experimental  MTF 


Figure  35.  2-D  indirect  photography  experimental  setup 


Four  objects  were  used  to  create  the  2-D  indirect  images:  (1)  two  white  1-cm 
squares  (each  square  subtending  24  mrad  with  respect  to  the  wall-to-object  distance) 
separated  by  1  cm  both  horizontally  and  vertically;  (2)  two  white  5 -mm  squares  (12 
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mrad)  separated  by  5  mm  both  horizontally  and  vertically;  (3)  two  white  2-mm  squares 
(4.8  mrad)  separated  by  2  mm  both  horizontally  and  vertically;  (4)  a  white  25-mm  (60 
mrad)  square  with  a  5 -mm  (12  mrad)  square  cut  from  the  center.  The  objects  were 
created  from  white  cardstock  on  a  flat  black  poster  board  background.  Figure  36  shows 
object  2  and  object  4.  The  semi-gloss  wall  and  imaged  reflectors  from  the  1-D 
experiment  were  used. 


(a)  (b) 

Figure  36.  2-D  (a)  object  2  and  (b)  object  4 


To  create  the  indirect  images  of  objects  1,  2  and  3,  2601  digital  images  in  a  51x51 
pattern  were  acquired  with  the  object  translated  horizontally  and/or  vertically  one  fifth  of 
the  square  size  between  each  data  image,  i.e.  2-mm  movement  for  object  1,  1-mm 
movement  for  object  2.  Three  indirect  images  were  created  of  object  4:  (1)  2601  digital 
images  in  a  51x51  pattern  with  the  object  translated  1-mm  horizontally  and/or  vertically 
between  each  data  image;  (2)  441  digital  photographs  in  a  21x21  pattern  with  2.5-mm 
movement  between  data  images;  and  (3)  121  digital  images  in  an  11x11  pattern  with  5- 
mm  movement  between  each  data  image,  which  corresponds  to  the  Nyquist  frequency  for 
object  4.  For  reference,  Figure  37  shows  an  overlay  of  the  reflected  laser  spot  and  object 
2,  and  Figure  38  shows  the  complete  2-D  laboratory  setup. 
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Figure  38.  2-D  laboratory  setup 
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As  with  the  1-D  dual  and  indirect  images,  to  improve  the  signal-to-noise  ratio,  the 
entire  digital  image  was  used  to  create  the  raw  indirect  image.  Initially,  the  average 
intensity  of  the  image's  pixels  was  used  to  form  the  intensity  on  the  imaging  system  for 
each  position,  x',  as  represented  by  Eq.  (83).  Figure  39  shows  the  raw  indirect  images  of 
objects  1,  2  and  3  and  the  indirect  images  of  object  4  created  with  2601,  441,  and  121 
images. 
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(d) 


(f) 


(e) 

Figure  39.  Unimproved  indirect  images  of  (a)  object  1,  (b)  object  2,  (c)  object  3,  and 
object  4's  (d)  51x51  image,  (e)  21x21  image  and  (f)  11x11  image 

Image  quality  was  then  improved  in  two  parts.  The  first  consisted  of  creating  a 
cumulative  histogram  of  each  recorded  image,  assigning  to  each  position,  x' ,  the  intensity 
corresponding  to  the  99.5th  percentile  of  the  cumulative  histogram,  i.e.  the  intensity  at 
which  99.5%  of  the  pixels  are  below  that  point  and  0.5%  are  above  that  point.  The 
subsequent  image  was  then  improved  using  MATLAB's  blind  deconvolution  command 
deconvblind.  The  initial  point  spread  function  for  the  deconvolution  was  a  block  of  ones 
one  pixel  less  than  the  image  size,  i.e.  50x50  for  the  51x51  image,  20x20  for  the  21x21 
image,  etc.  The  deconvolution  was  allowed  to  run  from  four  to  40  iterations  in  blocks  of 
4  iterations.  Figure  40  shows  the  averaged  51x51  indirect  image  of  object  four 
deconvolved  using  four,  eight,  16,  24,  32  and  40  iterations.  Figure  41  shows  the 
corresponding  cumulative  51x51  indirect  images  of  object  4. 
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(a)  (b)  (c) 


(d)  (e)  (f) 

Figure  40.  51x51  averaged  object  4  indirect  image  deconvolved  for  (a)  4  (b)  8  (c)  16 

(d)  24  (e)  32  and  (f)  40  iterations 

Following  each  block  of  four  iterations  in  the  deconvolution  process,  image 
quality  for  each  image  was  assessed  using  Eq.  (36).  Because  it  was  not  guaranteed  the 
object  would  be  in  the  exact  center  of  the  indirect  image,  every  possible  position  of  the 
object  was  evaluated  and  the  estimated  position  of  the  object  was  assigned  where  the 
image  quality  was  the  highest.  Figure  42  shows  the  image  quality  for  every  four 
iterations  of  the  deconvolution  process  for  the  51x51  cumulative  indirect  image  of  object 
4,  while  Figure  43  shows  the  highest  quality  image  produced  by  the  deconvolution 
process  and  the  overlay  of  the  ideal  image. 
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Figure  41.  51x51  cumulative  object  4  indirect  image  deconvolved  (a)  4  (b)  8  (c)  16 

(d)  24  (e)  32  and  (f)  40  iterations 

The  same  average  and  cumulative  techniques  and  deconvolution  processes  were 
accomplished  for  each  of  the  unimproved  indirect  images  seen  in  Figure  39.  The 
recovered  images  with  the  highest  image  quality  for  objects  1,  2  and  3  and  the  three 
indirect  images  of  object  4  are  shown  in  Figure  44. 
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Figure  42.  Image  quality  as  a  function  of  number  of  deconvolution  iterations  for 

51x51  indirect  image  of  object  4 


(a) 


(b) 


Figure  43.  (a)  Best  recovered  (deconvolved)  cumulative  image  of  object  4  and  (b) 
the  same  image  with  ideal  image  of  object  4  overlaid. 
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(C) 


O')  k)  (1) 

Figure  44.  Recovered  averaged  indirect  images  of  (a)  object  1,  (b)  object  2  and  (c) 
object  3.  Cumulative  indirect  images  of  (d)  object  1,  (e)  object  2  and  (f)  object  3. 
Averaged  (g)  51x51,  (h)  21x21  and  (i)  11x11  and  cumulative  (j)  51x51, 

(k)  21x21  and  (1)  11x11  indirect  images  of  object  4. 
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The  image  quality  analysis  of  objects  1,  2  and  3  can  be  converted  to  an  MTF  by 
considering  the  size  of  the  square  used  to  make  each  image  to  be  a  half  cycle.  Figure  45 
shows  the  MTF  for  unimproved  average  and  cumulative  indirect  images  and  the 
corresponding  deconvolved  indirect  images  of  objects  1,  2  and  3. 


0.05  o'l  0.15  0.2  0.25 

Cycles/mrad 


Figure  45.  MTF  of  unimproved  and  deconvolved  images  of  objects  1,  2  and  3. 

While  the  images  of  object  4  do  not  lend  themselves  to  a  traditional  MTF,  i.e. 
percent  modulation  or  image  quality  v.  cycles/(m)rad,  they  can  be  used  to  determine  the 
image  quality  v.  step  size  used  in  comparison  to  the  object  feature  size.  To  that  end,  the 
image  quality  of  the  three  indirect  images  and  the  associated  recovered  images  of  object  4 
were  plotted  as  function  of  ratio  of  the  step  size  to  feature  size,  i.e.  1-mm  step  size  and  5- 
mm  square  feature  on  object  4  yields  the  ratio  0.2.  Figure  46  shows  the  results  of  the 
analysis  for  object  4's  images. 

Indirect  images  of  object  2  were  accomplished  in  a  21x21  pattern  using  a  2.5-mm 
step  size  and  an  11x11  pattern  using  a  5-mm  step  size.  As  with  the  11x11  image  of 
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object  4,  the  1  lxl  1  image  of  object  2  represents  the  Nyquist  frequency  of  object  2.  The 
resultant  plot  of  image  quality  v.  the  ratio  of  step  size  per  feature  size  is  shown  in  Figure 
47. 


Step  per  Feature  Size 

Figure  46.  Object  4  image  quality  v.  step  size  results  for  51x51,  21x21  and  11x11 

images 

The  goal  of  the  2-D  setup  was  to  validate  that  indirect  photography  could  be 
expanded  from  the  1-D  sinusoidal  slides  with  well-defined  frequencies  to  more 
representative  2-D  objects.  Figure  44  shows  the  recovered  images  can  be  recognized; 
even  the  2-mm  (4.8  cycles/mrad)  squares  can  be  recognized  as  separate  squares  after  the 
deconvolution.  More  importantly,  Figure  45  shows  the  image  quality  can  be  improved 
by  a  blind  deconvolution. 
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Step  per  Feature  Size 

Figure  47.  Object  2  image  quality  v.  step  size  results  for  51x51,  21x21  and  11x11 

images 

While  the  experimental  setup  was  consistent  with  the  setup  used  in  the 
radiometric  model,  in  order  to  be  operationally  significant,  the  camera  and  the  laser  must 
be  co-located.  The  next  section  describes  the  results  of  co-locating  the  camera  and  the 
laser. 

Real-world  Setup 

Following  the  completion  of  the  2-D  experimental  setup,  the  camera  was  co¬ 
located  with  the  laser  creating  a  real-world  setup.  Figure  48  shows  the  experimental 
setup  and  Figure  49  shows  the  real-world  laboratory  setup.  While  the  separated  wall 
reflector  was  used,  it  was  placed  parallel  to  the  imaged  reflector  to  simulate  they  were 
part  of  the  same  wall. 
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Wall  Reflector  (w) 


Laser 


Figure  48.  Real  world  experimental  setup 


Initially,  object  2  from  the  2-D  experiment  described  in  the  previous  section  was 
used  to  test  the  real-world  configuration.  Figure  50  (a)  shows  the  dual  image  of  object  2, 
(b)  shows  the  indirect  image  and  (c)  shows  the  deconvolved  image.  The  indirect  image 
(b)  shows  a  banding  on  the  left  half  of  the  image  which  carries  over  to  the  deconvolved 
image.  It  was  determined  that  the  banding  was  caused  by  the  translation  stage.  As  the 
object  was  translated  to  the  left,  the  black  cardstock  covered  the  comer  of  the  translation 
stage  and  prevented  the  reflections  from  the  translation  stage  from  being  imaged  by  the 
camera  off  the  imaged  reflector.  Figure  5 1  shows  the  indirect  image  (a)  and  deconvolved 
image  (b)  with  the  block  in  place. 

Following  confirmation  the  2-D  setup  could  be  transitioned  to  a  real-world  setup, 
the  objects  were  changed  to  playing  cards.  Figure  52  shows  the  indirect  illumination  of 
the  five  of  clubs.  Raw  indirect  images  of  the  five  of  clubs  were  created  at  the  following 
resolutions:  (1)  96x63,  (2)  47x31,  (3)  23x15,  (4)  11x7  and  (5)  5x3.  Additionally,  a  5x3 
dual  image  was  created  (see  Figure  53). 
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Laser 
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Figure  49.  Real  world  laboratory  setup 


(a)  (b)  (c) 

Figure  50.  Object  2  real  world  setup  results 
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(a)  (b) 

Figure  51.  Real-world  (a)  raw  and  (b)  improved  indirect  images  of  object  1  with  the 

translation  stage  covered 


Figure  52.  Indirect  illumination  of  the  five  of  clubs 


At  this  point,  the  goal  of  the  experiment  was  to  identify  the  value  of  the  playing  card,  i.e. 
ace  through  king,  but  not  necessarily  the  suit.  To  that  end,  the  image  quality  of  the  1 1x7 
and  5x3  indirect  images  were  computed  as  though  they  were  the  ace  through  eight. 
(Note:  The  7  is  the  only  card  that  is  not  horizontally  symmetric;  therefore,  whether  the 
pip  is  in  the  upper  or  lower  position  must  be  tested  and  reported,  i.e.  7U  and  7L).  The 
resultant  image  qualities  are  reported  in  Table  6.  Based  on  the  results  presented  in  the 
previous  section,  the  image  quality  should  improve  as  the  step  size  to  feature  size 
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increases  (see  Figure  46  and  47).  Given  that  the  five  of  clubs  was  correctly  identified 
using  the  5x3  indirect  image,  the  next  set  of  experiments  will  be  used  to  try  to  correctly 
identify  the  ace  through  eight  using  5x3  indirect  images. 


(d)  (e)  (f) 

Figure  53.  5  of  clubs  indirect  images:  (a)  95x63,  (b)  47x31,  (c)  23x15,  (d)  11x7,  (e) 

5x3  and  (f)  5x3  dual  image 


Because  the  process  of  creating  indirect  images  is  both  data  and  time  intensive, 
using  the  lowest  resolution  possible  will  increase  the  operational  utility  of  an  indirect 
imaging  system.  To  that  end,  indirect  images  of  the  ace  of  clubs  through  the  eight  of 
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clubs  were  created  with  a  5x3  resolution  (see  Figure  54).  Following  the  creation  of  the 
indirect  images,  the  image  quality  of  each  image  was  calculated  for  each  of  the  possible 
cards.  The  highest  image  qualities  for  each  indirect  image  are  identified  in  red.  In  every 
case,  the  correct  card  was  identified. 


Table  6.  Five  of  clubs  indirect  image  quality 


Res 

Ace 

2 

3 

4 

Mask 

5 

6 

7L 

7U 

8 

5x3 

0.016 

0.001 

0.007 

0.026 

0.027 

0.015 

0.011 

0.012 

0.009 

11x7 

0.021 

0.012 

0.016 

0.032 

0.034 

0.016 

0.012 

0.013 

0.014 

Table  7.  Card  indirect  image  selection  criteria 


Card 

Ace 

2 

3 

4 

Mask 

5 

6 

7L 

7U 

8 

Ace 

0.032 

-0.005 

0.008 

-0.004 

0.004 

-0.002 

-0.002 

-0.002 

-0.002 

2 

-0.009 

0.032 

0.019 

0.003 

0.0002 

-0.003 

-0.004 

-0.003 

-0.005 

3 

0.022 

0.024 

0.026 

-0.001 

0.005 

-0.003 

-0.005 

-0.005 

-0.007 

4 

-0.015 

0.004 

-0.003 

0.032 

0.024 

0.018 

0.014 

0.015 

0.011 

5 

0.016 

0.001 

0.007 

0.026 

0.027 

0.015 

0.011 

0.012 

0.009 

6 

-0.007 

-0.004 

-0.005 

0.023 

0.018 

0.026 

0.021 

0.021 

0.017 

7 

-0.002 

-0.007 

-0.006 

0.015 

0.013 

0.018 

0.021 

0.013 

0.016 

8 

-0.005 

-0.009 

-0.008 

0.010 

0.008 

0.016 

0.014 

0.016 

0.018 

Following  the  identification  of  each  card,  the  image  quality  of  each  indirect  image 
was  improved  using  MATLAB's  blind  deconvolution  algorithm.  Figure  55  shows  the 
image  quality  improvement  per  iteration  for  the  two,  five  and  eight,  while  Table  8  lists 
the  image  quality  improvement  for  each  card  following  100  iterations  of  the 
deconvolution  algorithm  and  Figure  56  shows  the  improved  indirect  images. 
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Figure  54.  5x3  indirect  images  of  (a)  ace,  (b)  two,  (c)  three,  (d)  four,  (e)  five,  (f)  six, 

(g)  seven  and  (h)  eight  of  clubs 


- Two 


So  20  30  40  50  60  70  80  90  100 

Number  of  Deconvolution  Iterations 


Figure  55.  Image  quality  improvement  per  deconvolution  iteration 
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Table  8.  Indirect  image  improved  image  quality 


Card 

Card 

Ace 

2 

3 

4 

5 

6 

7L 

8 

Raw 

0.032 

0.032 

0.026 

0.032 

0.027 

0.026 

0.021 

0.018 

Improved 

0.902 

0.451 

0.317 

0.355 

0.329 

0.2641 

0.231 

0.2045 

Figure  56.  5x3  improved  indirect  images  for  the  (a)  ace,  (b)  two,  (c)  three,  (d)  four, 
(e)  five,  (f)  six,  (g)  seven  and  (h)  eight  of  clubs 
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Conclusion 


The  goal  of  the  simulation  and  experiments  was  to  confirm  the  radiometric  theory 
of  dual  and  indirect  photography  following  the  setups  shown  in  Figures  11  and  17.  In 
particular,  that  Eqs.  (70),  (83)  and  (96)  could  be  used  to  improve  the  quality  of  the 
images  as  defined  by  Eq.  (36). 

The  1-D  experiment  showed  indirect  photography  could  be  used  to  recover 
information  about  the  spatial-frequency  content  of  an  object  using  wall  and  imaged 
reflectors  with  different  BRDFs.  The  2-D  experiment  expanded  the  1-D  experiment  and 
produced  images  of  non-sinusoidal  objects.  The  real-world  experiment  confirmed 
indirect  photography  could  be  used  to  produce  images  in  an  operational  setup  by  imaging 
the  ace  through  eight  of  clubs  and  computationally  identifying  each  card  correctly. 

The  experiments  also  showed  the  image  quality  of  the  dual  or  indirect  images 
could  be  improved  though  a  deconvolution  process.  In  the  case  of  the  dual  experiment, 
the  small  angle-approximation  equation,  and  the  cross  section  of  the  laser  as  the  PSF  or  a 
blind  deconvolution  could  be  used  to  increase  the  image  quality.  In  the  case  of  indirect 
photography  a  blind  deconvolution  could  be  used  to  increase  the  image  quality. 

While  standard  deconvolution  techniques  did  increase  the  image  quality,  the 
symmetry  of  the  dual/indirect  photography  process  may  allow  for  further  improvement  of 
the  image  quality.  The  next  chapter  lays  the  foundation  for  that  further  improvement. 


81 


IV.  Matrix  Formulation 


While  standard  deconvolution  techniques  can  be  used  to  improve  the  image 
quality  of  the  indirect  images,  the  formation  of  the  indirect  images  creates  symmetries 
which  may  offer  the  opportunity  to  improve  the  deconvolution  process.  By  modeling  the 
indirect  image  equation  using  matrices,  some  of  these  symmetries  have  been  revealed  and 
research  is  continuing  in  the  effort  to  improve  the  image  quality. 


Matrix  theory  of  indirect  photography 

To  create  the  matrix  representation  of  indirect  photography,  each  component  of 
the  indirect  photography  equation,  Eq.  (78),  is  represented  by  a  matrix  resulting  in  the 
following  equation: 


ID)  =  (IB  ©  ¥im)  •  (A  ©  ¥ph)  -PO((GO  Fw)  •  Ew)  (98) 


where 


ID)  is  a  y  x  n  matrix  representing  the  data,  y  is  the  number  of  pixels  in  the  camera 
and  n  is  the  number  of  data  images. 

Ew  is  a  w  x  n  matrix  representing  the  irradiance  on  the  wall  reflector  and  w  is  the 
number  of  individual  points  on  the  wall. 

¥w  is  a  x  x  w  matrix  representing  the  BRDF  of  the  wall  from  every  point  on  the 
wall  to  every  point  in  the  fixed  object  frame  of  reference  and  x  is  the  total 
number  of  points  in  the  fixed  object  frame  of  reference. 

<G  is  a  x  x  w  matrix  representing  the  geometry  terms  (y)  from  every  point  on  the 
wall  to  every  point  in  the  fixed  object  frame  of  reference. 

P  ,  the  object  reflectance  matrix,  is  a  x  xn  matrix  representing  the  position  of  the 
object  of  interest  in  the  fixed  object  frame  of  reference  for  every  data  image. 

¥ph  is  a  y  x  x  matrix  representing  the  BRDF  from  the  fixed  object  frame  of 
reference  to  the  imaged  reflector  and  y  is  the  number  pixels  in  the  camera.3 

A  is  a  y  x  x  matrix  representing  the  geometry  terms  (a)  from  every  point  in  the 


3  The  points  on  the  imaged  reflector  correspond  to  the  projection  of  the  camera’s  pixels  onto  the  imaged 
reflector. 


82 


fixed  object  frame  of  reference  to  every  point  on  the  imaged  reflector. 

Wim  is  a  y  x  y  matrix  representing  the  BRDF  from  the  imaged  reflector  to  the 
lens  of  the  imaging  system  which  is  subsequently  focused  on  the  camera's 
individual  pixels.4 

IB  is  a  y  x  y  matrix  representing  the  geometry  terms  (/?)  from  every  point  in  the 
fixed  object  frame  of  reference  to  the  pixels  in  the  camera. 

•  represents  standard  matrix  multiplication. 

O  represents  the  Hadamard  product. 

With  the  goal  of  recovering  the  image  represented  by  any  column  of  the  object 
reflectance  matrix,  P,  the  matrices  before  and  after  the  reflectance  matrix  can  be 
evaluated  to  form  the  equation: 

®yxn  =  Qyxx  ‘  ^xxn  O  H^exn  (9^) 

where: 

Qyxx  =  (®  O  IFijn)yxy  *  O  I ^ph)yxx  (100) 

(GO 

F wall)xxw  *  ®wa^WXn  (101) 

To  identify  the  symmetry  created  by  the  indirect  image  process,  the  structure  of 
each  matrix  of  Eq.  (99)  will  be  evaluated.  Based  on  the  unknown  BRDFs  of  both  the 
imaged  reflector,  ¥im,  and  the  phase  function  of  the  object,  Wph,  as  well  as  the  unknown 
geometry  of  the  setup  between  the  object,  the  imaged  reflector  and  the  lens  of  the 
imaging  system,  nothing  can  be  definitely  stated  about  the  structure  of  the  Q  matrix 
without  a  priori  knowledge.  Therefore,  the  Q  matrix,  in  its  most  general  form,  is 
represented  by  a  y  x  x  matrix  of  unknown  elements. 


4  With  an  ideal  imaging  system,  Wim  would  be  a  diagonal  matrix,  i.e.  each  pixel  is  perfectly  focused. 
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'yxx  — 


(102) 


<7 11  <7ix' 

_9yl  9yx_ 

Since  both  the  dual  and  indirect  photography  algorithms  require  the  same  set  of  pixels 
from  each  data  image  to  be  used  to  create  the  image,  and  in  keeping  with  the  experiments 
detailed  in  Chapter  IV,  where  the  entire  image  is  used  to  form  the  dual/indirect  image,  the 
(Q>  matrix  can  be  represented  by  a  row  vector: 

Qzx*  =  [9 1’  -  ,qn]  (103) 

The  structure  of  the  object  reflectance  matrix,  P,  and  the  irradiance  matrix,  R,  is 
determined  by  the  number  of  data  images  taken,  n,  the  distance  the  object  of  interest  is 
translated  between  each  data  image  in  comparison  to  the  size  of  the  object  of  interest,  and 
the  pattern  in  which  the  object  is  translated.  To  demonstrate  this  concept,  an  object  with 
four  distinctive  points  in  a  two-by-two  square  pattern  will  be  used.  The  first  data  image 
is  acquired  for  an  indirect  image  as  described  in  the  real-world  section  of  Chapter  IV. 
The  object  is  then  translated  vertically  downward  a  distance  equal  to  one  half  the  vertical 
length  of  the  object  and  a  second  data  image  is  acquired.  For  the  second  data  image,  the 
irradiance  on  the  upper  left  quadrant  of  the  object  will  be  the  same  as  the  irradiance  on 
the  lower  left  quadrant  of  the  object  in  the  first  image.  Likewise,  the  irradiance  on  the 
upper  right  quadrant  of  the  object  in  the  second  data  image  is  the  same  as  the  irradiance 
on  the  lower  right  quadrant  in  the  first  data  image.  Figure  57  illustrates  this  symmetry 
where  Figure  57  (a)  is  the  position  of  the  object  when  the  first  data  image  is  acquired  in 
comparison  to  the  irradiance  in  the  fixed  object  frame  of  reference.  Figure  57  (b)  shows 
the  position  of  the  object  after  the  translation  vertically  downward.  In  both  cases,  the 
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blue  (upper  left)  designators  represent  the  irradiance  in  the  fixed  object  frame  of 
reference  and  the  red  (lower  right)  designators  represent  the  reflectance  of  that  quadrant 
of  the  object. 


Figure  57.  Irradiance  of  the  object  given  a  2x2  pattern  for  the  (a)  first,  (b)  second, 

(c)  third  and  (d)  fourth  data  images. 

Following  the  second  data  image,  the  object  is  translated  vertically  up,  to  the 
original  vertical  position  and  then  horizontally  to  the  right,  a  distance  equal  to  one  half 
the  horizontal  width  of  the  object  resulting  in  the  configuration  shown  in  Figure  57(c).  In 
this  position,  the  irradiance  on  the  upper  and  lower  left  quadrants  of  the  object  are  the 
same  as  the  irradiance  on  the  upper  and  lower  right  quadrants  of  the  object  in  the  first 
data  image.  Figure  57(d)  shows  the  position  of  the  object  in  the  fixed  object  frame  of 
reference  after  it  has  been  translated  vertically  downward  from  the  position  of  the  object 
in  the  third  data  image.  The  relationship  between  the  irradiance  in  the  fixed  object  frame 
of  reference  and  the  reflectance  of  the  object  between  the  third  and  fourth  data  images  is 
the  same  as  the  relationship  previously  described  between  the  first  and  second  data 
images. 

Given  the  fixed  geometry  between  the  wall  reflector  and  the  fixed  object  frame  of 
reference,  the  irradiance  in  the  fixed  object  frame  of  reference  is  unchanged  from  one 
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data  image  to  the  next.  Given  the  scenario  described  in  Figure  57,  the  irradiance  matrix, 


E,  can  be  represented  by: 

rr±  rx  ry  rr i 
r2  r2  r2  r2 

r3  r3  r3  r3 

r4  r4  r4  r4 

Exxn  =  r5  r5  r5  r5  (104) 


where  the  individual  elements,  r1(  r2,  etc,,  represent  the  irradiance  incident  on  the  specific 
area  in  the  fixed  object  frame  of  reference.  In  a  more  general  form,  the  irradiance  matrix, 
E,  can  be  written  as  the  Kronecker  product  of  two  vectors  R  and  e„,  where  R  is  the 
column  vector  that  results  from  applying  the  Vec  operator  on  the  matrix  describing  the 
irradiance  in  the  fixed  object  frame  of  reference  and  the  vector  en  is  a  column  vector  of 
n  ones  where  n  is  the  number  of  data  images,  (see  Eq.  (105)) 

E  =  R  ®  e„  (105) 


In  the  case  of  the  scenario  described  by  Figure  57,  R  and  en  are  represented  by: 

ri 

r2 

*3 
r4 

r5 
r6 

r7 

r8 

r9 

The  object  reflectance  matrix,  P,  is  an  x  xn  matrix  where,  for  each  column  of  the 
matrix,  the  Vec  operator  has  been  applied  to  the  matrix  describing  the  object  in  the  fixed 
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object  frame  of  reference  corresponding  to  the  data  image.  The  P  matrix  for  the  scenario 


described  in  Figure  57  is: 


Pi  0  0  0' 

P2  Pi  0  0 

0  p2  0  0 

p3  0  pi  0 

'xxn  =  P4  P3  P2  Pi 

o  p4  0  p2 

0  0  P3  0 

00  p4  Ps 
-  0  0  o  p4_ 


(107) 


where  the  rows  of  the  matrix  correspond  to  the  distinctive  points  in  the  fixed  object  frame 
of  reference  and  the  columns  correspond  to  the  different  data  images.  Substituting  Eqs. 
(103),  (104)  and  (107)  into  Eq.  (99),  the  equation  for  the  data  matrix  becomes: 


(\Pi 

P2 

0 

Ps 

Dl xn  =  [9i  -  <79]  •  p4 

0 

0 

\lo 


0  0  0  ]  -ry  ly  ry 

Pi  0  0  r2  r2  r2 

P2  0  0  r3  r3  r3 

0  Pi  0  r4  r4  r4 

Ps  P2  Pi  O  rs  r5  r5 

p4  0  p2  r6  r6  r6 
0  p3  0  r7  r7  r7 

0  p4  Ps  rs  rs  rs 

0  0  p4J  Lr9  r9  r9 


(108) 


As  previously  stated,  the  transport  of  light  through  the  system  described  above 
and  used  to  create  dual/indirect  images  is  linear.  Given  that  the  E  and  Q  matrices  are 
defined  by  the  geometry  of  the  setup  and  the  irradiance  of  the  laser  spot,  both  constant 
throughout  the  creation  of  the  indirect  image,  by  creating  two  basis  sets,  23 P  and  to 
describe  the  reflectance  matrix,  P,  and  the  data  matrix,  D,  respectively,  and  defined  as: 

=  {Y1(Y2,Y3,Y4} 
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(109) 

(110) 


where 
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and 


1 

-1- 

1  1 

-0- 

1  1 

o 

II 

0 

0 

II 

1 

0 

II 

0 

0 

1 

-0- 

1  1 

-0- 

1  1 

-1- 

(112) 


then  Eq.  (108)  can  be  modeled  as  a  linear  transform  from  the  object  reflectance  basis  to 
the  data  basis  formally  defined  as: 

£(P,  D)  =  {X:  P  ->  D]  (113) 


where  the  linear  transform,  %,  is  formed  by  applying  the  Vec  operator  to  the  data  matrix, 
D,  created  when  the  respective  object  reflectance  basis  sets  are  evaluated.  For  the 
scenario  described  in  Figure  57 

%  =  [dm,  dm,  dm,  ©M]  0 14) 

Evaluating  Eq.  (114),  the  linear  transform  becomes: 


qiri 

q2r2 

q4r4 

qs^' 

qzr2 

qsr3 

q5r5 

qe^e 

q4r4 

qsr5 

q7r7 

qs^s 

.qsr5 

qer6 

qsr8 

q9r9_ 

(115) 


which  is  a  block-Hankel,  Hankel-block  matrix  and  can  be  simplified  to: 
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"a  b  c  d" 

~  _  b  e  d  f 

c  d  g  h 

,d  f  h  i. 

where 


a  =  q^,  b  =  q2r2, ,  i  =  q9r9 


(116) 


(117) 


Using  Eq.  (116),  an  indirect  image  can  be  modeled  as  the  linear  transform  operating  on  a 
column  vector  representing  the  object  to  produce  the  recorded  data: 
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In  keeping  with  the  theory  of  indirect  photography  described  in  Chapter  III,  if  the 
BRDF  and  geometry  of  the  setup  are  known  a  priori,  then  the  transform  matrix  will  be 
known.  The  transform  matrix  can  then  be  inverted  and  when  both  sides  of  Eq.  (118)  are 
multiplied  from  the  left  by  the  inverted  transform  matrix,  XT1,  the  reflectance  matrix,  IP, 
can  be  solved  for,  which,  in-tum,  allows  for  the  reconstruction  of  the  object  of  interest. 

Given  that  the  BRDFs  and  geometry  of  the  setup  will  likely  not  be  available  in  an 
operational  environment,  Eq.  (118)  must  be  solved  without  explicitly  knowing  the 
transform  matrix,  X.  Eq.  (118)  is  underspecified,  with  four  equations  and  13  unknowns. 
Therefore,  solving  the  system  of  equations  directly  will  not  be  possible.  However,  it  may 
be  possible  to  solve  Eq.  (118)  by  posing  it  as  an  optimization  problem  and  finding  the 
optimum  transform  matrix  and  optimum  image  vector  given  the  constraint  that  the 
multiplication  of  the  two  matrices  results  in  the  data  matrix. 
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To  that  end,  the  object  has  been  simplified  to  a  2x1  matrix,  resulting  in  a  2x2 
transform  matrix  and  a  2x1  data  matrix.  This  results  in  two  equations  and  five 
unknowns.  While  research  to  date  has  not  yielded  a  solution,  progress  has  been  made  by 
optimizing  the  solution  such  that  total  energy  of  the  system  is  minimized  while  still 
satisfying  Eq.  (118). 

Conclusion 

At  the  time  of  this  writing,  the  matrix  formulation  has  yet  to  yield  the  desired 
results.  While  more  research  is  required,  I  believe  the  best  path  forward  is  to  divide  the 
data  images  into  equal  sections,  i.e.  halves,  quadrants,  etc,,  and  using  each  section  of  the 
recorded  data  to  form  an  indirect  image.  While  each  of  these  indirect  images  will  have  a 
unique  linear  transform,  the  object  remains  the  same  for  all  of  the  indirect  images. 
Therefore,  solving  Eq.  (118)  simultaneously  for  all  of  the  indirect  images  may  yield  a 
unique  solution  at  the  intersection  of  the  sets  of  solutions  formed  by  the  individual 
images. 


90 


VI.  Conclusion 


While  techniques  to  image  objects  through  triple  layered  jungle  canopies  and 
camouflage  netting  have  been  previously  developed,  the  theory  described  in  Chapter  III 
and  the  experimental  proof  of  concept  described  in  Chapter  IV  are  the  first  to  allow 
images  to  be  created  either  around  comers  or  of  objects  under  solid  shelters. 

The  concept  of  dual  photography,  originally  designed  to  aid  in  the  creation  of 
computer  generated  graphics,  was  radiometrically  modeled  and  simplified  to  Eq.  (70), 
repeated  here  as  Eq.  (119),  which  revealed  the  dual  image  was  a  convolution  of  the  object 
of  interest  and  a  kernel  comprised  of  the  reflection  characteristics  of  the  object  and  non- 
specular  reflectors  as  well  as  the  geometry  of  the  set  up  used  to  create  the  dual  image. 

®pixi(%  >  @obj') 


I  r i(x' -  x" ,  Gwaii)  p(x")  dx" 
Jlas 


(119) 


The  dual  photography  radiometric  equation  was  further  simplified  by  using  the  small- 
angle  approximation  to  Eq.  (95),  repeated  here  as  Eq.  (120),  which  revealed  the  dual 
image  could  be  approximated  as  a  convolution  of  the  object  of  interest  and  the 
illumination  source. 


O 


piXi 


(x)  =  Kj  I  Eobj(x'  —  x  )p(x  )dx 
Jlas 


(120) 


These  equations  showed  the  image  quality  of  the  dual  images  could  be  improved 
through  a  deconvolution  process.  In  the  case  of  the  small-angle  approximation 
Eq.  (120),  the  irradiance  on  the  object  of  interest  could  be  used  to  approximate  the  point 
spread  function  (PSF)  of  the  system  and  image  quality  could  be  improved  by  standard 
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deconvolution  techniques.  In  the  case  of  the  standard  dual  photography  equation 
Eq.  (119),  two  processes  could  be  used  to  improve  the  image  quality.  If  the  exact  BRDFs 
of  the  reflecting  surfaces  and  geometry  of  the  setup  are  known,  standard  deconvolution 
techniques  can  be  used  to  improve  the  image  quality.  The  second  process  used  a  blind 
deconvolution  to  improve  the  image  quality.  This  technique  has  the  advantage  of  not 
requiring  explicit  knowledge  of  the  geometry  and  BRDFs  of  the  reflecting  surfaces,  a 
constraint  likely  to  be  encountered  in  an  operational  situation. 

Following  the  development  of  the  dual  photography  radiometric  equations,  the 
irradiance  on  the  object  of  interest  was  changed  from  the  laser  spot  of  dual  photography 
to  a  reflection  from  a  non-specular  surface,  resulting  in  Eq.  (83)  repeated  here: 

As  with  the  dual  photography  equation,  the  indirect  photography  equation  implied  two 
important  concepts:  (1)  an  indirect  image  could  be  created  if  an  individual  pixel,  any 
group  of  pixels  or  the  entire  image  was  used  to  create  the  image,  as  long  as  the  same  sets 
of  pixels  were  used  from  every  data  image,  and  (2)  the  image  quality  of  the  resultant 
indirect  image  could  be  improved  through  a  (blind)  deconvolution  technique. 

Following  the  development  of  the  theoretical  equations,  the  dual  photography 
equations  were  experimentally  validated  using  six  sinusoidal  slides  as  objects.  Figure  28, 
repeated  here  as  Figure  58  (a),  shows  the  unimproved  dual  image  created  using  a  1.0 
cycles/mm  (0.20  cycles/  mrad)  slide  and  image  quality  improvements  made  using  Lucy- 
Richardson  and  blind  deconvolution  algorithms.  Figure  32,  repeated  here  as  Figure  58 
(b),  is  the  modulation  transfer  function  (MTF)  created  by  the  dual/indirect  images  of  all 


I  Tj(x'  -  x",  ewall)  p{x ")  dx" 

* obj 


(121) 
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six  slides  and  shows  the  improvement  in  image  quality  achieved  by  the  deconvolution 


process. 


45  47  49  51  53 

milli-meters 


(a)  (b) 

Figure  58.  Summary  of  1-D  experiment  (a)  Dual  1.0  cycles/mm  image  and 
deconvolutions  and  (b)  1-D  semi-gloss  reflector  MTF 


The  experiment  was  expanded  to  2-D,  resulting  in  indirect  images  being  created 
of  simple  geometric  objects.  The  resultant  image  quality  was  evaluated  for  the  raw  and 
improved  (deconvolved)  images.  Representative  images  showing  the  raw  and 
deconvolved  2-D  indirect  images  and  the  resultant  MTF  originally  shown  in  Figure  43 
and  45  have  been  reprinted  here  as  Figure  59. 

Finally,  the  experiment  was  reconfigured,  co-locating  the  camera  and  the  laser  in 
an  operationally  representative  configuration.  Indirect  images  of  eight  playing  cards 
were  created  and  evaluated  against  the  ideal  image  of  all  eight  playing  cards.  The  image 
quality  analysis  resulted  in  all  eight  playing  cards  being  properly  identified  at  a  resolution 
of  5x3  pixels.  The  indirect  images  of  the  5  of  clubs  at  various  resolutions  are  shown  in 
Figure  53  for  comparison  and  the  47x3 1  resolution  is  repeated  here  as  Figure  60. 
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Figure  59.  Indirect  image  of  a  square  annulus  (a)  unimproved  and  (b)  deconvolved 
ideal  image  of  the  annulus  overlaid  and(c)  MTF  created  from  offset  squares. 
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Figure  60.  47x31  indirect  image  of  the  5  of  clubs. 
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The  results  of  this  research  have  been  presented  at  the  SPIE  conference  on 
Reflections,  Scattering  and  Diffraction  from  Surfaces  II.  It  has  also  been  submitted  to 
Optics  Express  and  is  currently  under  peer  review  for  publication.  Additionally,  a  patent 
application  has  been  filed  with  Air  Force  Material  Command  for  the  concept  of  indirect 
photography  and  is  under  review  by  intellectual  property  lawyers. 

While  the  research  presented  in  this  document  achieved  the  desired  goal  of 
developing  the  theory  of  indirect  photography  and  experimentally  proving  that  theory, 
additional  improvements  to  the  theory  lend  themselves  to  follow-on  research.  Some  of 
these  possible  improvements  are:  (1)  the  continuation  of  the  concept  described  in 
Chapter  V,  where  the  improvement  to  the  deconvolution  process  is  sought  by  taking 
advantage  of  the  known  symmetries  in  the  matrices  that  emerge  from  the  radiometric 
equations  involved  in  the  creation  of  an  indirect  image.  (2)  The  second  area  of  possible 
continued  research  is  to  polarimetrically  model  the  dual  and  indirect  setups  in  an  effort  to 
take  advantage  of  reflectance  differences  between  polarizations  in  the  creation  of  the 
indirect  image.  (3)  Application  of  advanced  signal/image  processing  beyond  the 
deconvolution  process.  (4)  Research  into  removing  the  limitations/assumptions  of  the 
theory  described  in  Chapter  III,  i.e.  3-D  objects  and  objects  with  object  with  varying 
phase  functions  as  well  as  reflectance. 

The  technique  of  indirect  photography  described  in  this  document  is  still  in  the 
early  stages  of  development  and  requires  additional  research  before  an  operational 
prototype  can  be  fielded.  That  said,  I  believe  this  document  lays  the  foundation  for  that 
research  and  has  the  potential,  when  fully  developed,  to  aid  the  military  and  intelligence 
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communities  in  their  ability  to  identify  and  classify  items  of  interest  in  situations 
currently  not  possible. 
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Appendix  A.  Play  card  BRDFs 


To  validate  the  assumption  that  the  BRDF  of  playing  cards  can  be  decomposed 
into  a  reflectance,  p,  and  a  phase  function  that  controls  the  underlying  angular  shape  of 
the  BRDF  as  described  in  Eq.  (50),  AFIT’s  CASI  was  used  to  measure  the  BRDF  of  the 
white,  black  and  red  portions  of  a  standard  playing  cards.  Figure  61  shows  the  measured 
BRDFs  resulting  from  633nm  HeNe  laser  and  a  45  degree  incident  angle.  In  this  graph, 
‘X’  axis  is  the  angular  difference  from  the  specular  reflection  i.e.  0  on  the  ‘X’  axis 
represents  specular  reflection,  positive  angles  are  away  from  the  surface  normal  and 
negative  numbers  are  toward  the  surface  normal  and  the  incident  irradiance.  The  gap  in 
the  data  at  -90  degrees  comes  from  the  sensor  blocking  the  incident  irradiance. 
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Appendix  B.  Indirect  Photography  Simulation  MATLAB  Code 


The  following  MATLAB  code  was  used  to  simulate  the  dual  and  indirect 
photography  experiments.  A  wire-frame  model  of  the  fixed  object  frame  of  reference 
(25x25  reflector  points)  and  the  wall,  for  indirect  photography,  and  imaged  reflectors 
(3x3  and  21x21  reflector  points  respectively)  were  created.  Arrays  are  then  created  for 
all  incident  and  reflected  angles  and  distances  from  each  point  on  the  adjacent  reflectors. 
The  glint  angle  and  subsequent  BRDF  using  Eq.  (27)  on  page  18  are  also  created  and 
stored  in  arrays.  The  object  was  then  placed  in  the  upper  left  comer  of  the  fixed  object 
frame  of  reference  and  beginning  from  the  wall  reflector,  every  possible  path  to  the 
imaging  system  was  evaluated,  i.e.  from  each  point  on  the  wall  reflector  to  every  point  on 
the  object  reflector  to  every  point  on  the  imaged  reflector,  etc.  The  total  flux  impinging 
on  the  lens  of  the  imaging  system  was  summed  to  simulate  the  entire  image  was  being 
used  to  create  the  dual  or  indirect  image.  The  object  was  then  translated  horizontally  and 
vertically  through  each  of  the  possible  positions  and  the  process  repeated.  The  resultant 
images  are  then  improved  using  both  Lucy-Richardson  and  blind  deconvolution 
algorithms.  The  raw  and  improved  images  are  shown  in  Figure  20  and  22  on  pages  48 
and  50  respectively. 
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% . . - . - 

%This  code  is  a  dual  and  indirect  simulation  designed  to  test  the  basic 
%equations  and  assumptions  made  in  the  mathematical  model. 

% - 


% . — - 

%Required  Input 
% . .  . 

laser  =  [0,0,-l]; 

WallNormal  =  [0,1,-1]; 
WallCenterPoint  =  [0,0,0]; 
xNumWallPoints  =3; 
yNumWallPoints  =3; 
xWallLength  =  .08; 
yWallLength  =  .08; 

xNumWalllmPoints  =2 1 ; 
yNumWalllmPoints  =2 1 ; 
xWalllmLength  =  3; 
yWalllmLength  =  5; 

ObjNormal  =  [0,-l,0]; 
ObjCenterPoint  =  [0,5,0]; 
xNumObjPoints  = 25 ; 
zNumObjPoints  = 25 ; 
xObj  Length  =  1; 
zObj  Length  =  1; 

LensNormal  =  [0,0,1]; 
LensCenterPoint  =  [0,0,- 10]; 
xNumLensPoints  =3; 
yNumLensPoints  =3; 
xLensLength  =  1 ; 
yLensLength  =  1 ; 

SigmaWall  =  .5*  pi/180; 
SigmaObj  =  1.5  *  pi/180; 

xSlideWidth  =  1; 
zSlideWidth  =  1; 
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% . — - - -  - - 

%  Basic  Calculations 

WallNormalHat  =  WallNormal/norm(WallNormal); 
ObjNormalHat  =  ObjNormal/norm(Obj  Normal); 
LensNormalHat  =  LensNormal/norm(LensNormal); 

TotWallPoints  =  xNumWahPoints*yNumWallPoints; 
TotWalllmPoints  =  xNurnWahImPoints*yNumWahImPoints; 
TotObjPoints  =  xN  umObj  P  o  ints  *  zN  umObj  P  o  ints ; 
TotLensPoints  =  xNumLensPoints*yNumLensPoints; 

xSlideHalfWidth  =  xSlide  Width/2; 
zSlideHalfWidth  =  zSlide  Width/2; 


% - 

%  This  section  creates  the  Reflector  Points 
% - 


% . - . . 

%  This  section  creates  the  Wall  Reflector  Points 

xWallStart  =  WallCenterPoint(l,l)-  (xWallLength/2); 

yWallStart  =  WallCenterPoint(l,2)-  (yWahLength/2)*WahNormalHat(l,2); 

xWallStep  =  (xWallLength  /(xNumWallPoints-1)) ; 

yWallStep  «  (yWallLength  /(yNumWallPoints-l))*WallNormalHat(l,2); 

Wall  Array  =  zeros(xNumWahPoints,yNumWallPoints,3); 

for  i  =  1  :xNumWallPoints 

WallArray(i, : ,  1  )=  xWallStart+xWallStep*(i- 1 ); 
for  j  =  l:yNumWallPoints 

WallArray(:,j,2)=  yWahStart+yWallStep*(j-l); 

WallArray(:,j,3)=  (WahNormalHat(l,2)*(WallCenterPoint(l,2)-... 

(y W  allStart+y W  allStep  *(j  - 1  ))))/WallNormalHat(  1,3); 
end 
end 

FlatWallArray  =  squeeze(reshape(WallArray,l, TotWallPoints, 3)); 


% . . - . — . — .  .  . 

%  This  section  creates  the  Imaged  Wall  Reflector  Points 
xWalllmStart  =  WallCenterPoint(l,l)-  (xWallImLength/2); 
yWalllmStart  =  WallCenterPoint(l,2)-  (yWallImLength/2)*WallNormalHat(l,2); 

xWalllmStep  =  (xWalllmLength  /(xNumWalllmPoints-l))  ; 

yWalllmStep  =  (yWalllmLength  /(yNumWallImPoints-l))*WallNormalHat(l,2); 

Wall  I  m  Array  =  zeros(xNumWallImPoints,yNumWallImPoints,3); 

for  i  =  1  :xNumWahImPoints 

W  allImArray(i, : ,  1  )=  x  W  alllmS  tart+xW  alllmS  tep  *  (i  - 1 ) ; 
for  j  =  l:yNumWallImPoints 
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WallImArray(:,j  ,2)=  yWallImStart+yWallImStep*(j  - 1 ); 
WallImArray(:  ,j  ,3)=  (WallNormalHat(  1 ,2)  *(WallCenterPoint(  1,2)-... 
(yWallImStart+yWallImStep*(j  - 1  ))))/WallNormalHat(  1,3); 
end 
end 

FlatWalllmArray  =  squeeze(reshape(WallImArray,  1  ,TotWallImPoints,3)); 


% - 

%  This  section  creates  the  Object  Reflector  Points 
xObjStart  =  ObjCenterPoint(l,l)-  xObjLength/2; 
zObjStart  =  ObjCenterPoint(l,3)-  zObjLength/2; 

xObjStep  =  xObjLength  /(xNumObjPoints-1); 
zObjStep  =  zObjLength  /(zNumObjPoints-1); 

Obj  Array  =  zeros(xNumObj  Points, zNumObjPoints,3); 

for  i  =  1  :xNumObjPoints 
for  j  =  l:zNumObj  Points 

ObjArray(i,:,l)  =  xObjStart  +  xObjStep*(i-l); 

Obj  Array  (:,j,  3)  =  zObjStart  +  zObjStep*(j-l); 
end 
end 

ObjArray(:,:,2)  =  ObjCenterPoint(l,2); 

FlatObj  Array  =  squeeze(reshape(ObjArray,l,TotObjPoints,3)); 


% - - - — . — . - . - .  - . 

%  This  section  creates  the  Lens  Points 
xLensStart  =  LensCenterPoint(l,l)-  (xLensLength/2); 
yLensStart  =  LensCenterPoint(l,2)-  (yLensLength/2); 

xLensStep  =  (xLensLength  /(xNumLensPoints-1)); 
yLensStep  =  (yLensLength  /(yNumLensPoints-1)); 

Lens  Array  =  zeros(xNumLensPoints,yNumLensPoints,3); 

for  i  =  1  :xNumLensPoints 
for  j  =  1  :yNumLensPoints 

Lens  Array(i, : ,  1  )=xLensStart+xLensStep  *(i- 1 ) ; 

Lens  Array(:  ,j  ,2)=yLensStart+yLensStep*(j  - 1 ); 
end 
end 

LensArray(:,:,3)  =  LensCenterPoint(l,3); 

FlatLensArray  =  squeeze(reshape(LensArray,l,TotLensPoints,3)); 
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clear  (’WallCenterPoint’,’xWallLength’,’yWallLength’,... 
'xW  alllmLength’ ,  ’y  W  alllmLength’ , . . . 

’Obj  CenterPointVxObj  Length ’,’zObj  Length’, • 
’LensCenterPointVxLensLengthVyLensLength’) 

clear  (’xLensStartVxLensStepVyLensStartVyLensStep’,... 
’xWallImStart’,’xWallImStep’,’xWallStart’,’xWallStep’,... 
'yWallImStart’,'yWallImStep’,’yWallStart','yWallStep',... 
’xObj  Start’,’xObj  Step’,’zObj  Start’,’zObj  Step’) 

clear('Lens  Array’,’ Wall  Array’, ’WalllmArray’) 

clear(’WallNormar,’ObjNormar,’LensNormar) 


% - 

%  This  section  creates  the  Distance  Arrays  for  calculations 
% - 


% - - - - - - - - 

%  Distance  from  every  point  on  the  Wall  to  every  point  on  the  Object 
DistToObj  =  zeros(TotWallPoints,TotObjPoints); 

for  i=l:TotWallPoints 
for  j = 1 :  T  otObj  P  oints 

DistT  oObj  (i,j  )=  norm(FlatW  all  Array(i, : )  -FlatObj  Array(j , : )) ; 
end 
end 


%  Distance  from  every  point  on  the  Obj  to  every  point  on  the  Imaged  Wall 

DistToWalllm  =  zeros(TotObjPoints,TotWallImPoints); 

for  i=l:TotObjPoints 
for  j=l  :TotWallImPoints 

DistToWallIm(i,j)=  norm(FlatObjArray(i,:)-FlatWallImArray(j,:)); 
end 
end 


% - 

%  Distance  from  every  point  on  the  Imaged  Wall  to  every  point  on  the  Lens 

DistToLens  =  zeros(TotWallPoints,TotLensPoints); 

for  i=l:TotWallImPoints 
for  j=l  :TotLensPoints 

DistT  oLens(i,j  )=  norm(F  latW  allImArray(i, : )  -FlatLens  Array(j , : )) ; 
end 
end 
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% . . — . — . — . — . - 

%  This  section  creates  the  Theta  Arrays  for  calculations 
% . . - . . . - . ----- 


% . --- . - . --- . - 

%Cos  of  the  angle  between  the  Wall  Normal  and  every  vector  between  the  Wall 
%points  and  the  Object  Points 

CosThetaWallPrime  =  zeros(TotWallPoints,TotObjPoints); 

for  i=l:TotWallPoints 
for  j = 1 :  T  otObj  P  oints 

RefVect  =  FlatObjArray(j,:)-FlatWallArray(i,:); 

RefVectHat  =  RefVect/norm(RefVect); 

CosThetaWallPrime(i,j)=  dot(WallNormalHat,RefV  ectHat); 
end 
end 


% . - . 

%Cos  of  the  angle  between  the  Obj  Normal  and  every  vector  between  the  Wall 
%points  and  the  Object  Points 

CosThetaObj  =  zeros(TotWallPoints,TotObjPoints); 

for  i=l:TotWallPoints 
for  j=l:TotObjPoints 

IncVect  =  FlatWallArray(i,:)-FlatObjArray(j,:); 

IncVectHat  =  IncVect/norm(IncVect); 

CosThetaObj  (i,j)=  dot(ObjNormalHat, IncVectHat); 
end 
end 


% - 

%Cos  of  the  angle  between  the  Obj  Normal  and  every  vector  between  the 
%Imaged  Wall  points  and  the  Object  Points 

CosThetaObjPrime  =  zeros(TotObjPoints,TotWallImPoints); 

for  i=l:TotObjPoints 
for  j=l  :TotWallImPoints 

RefVect  =  FlatWallImArray(j,:)-FlatObjArray(i,:); 

RefVectHat  =  RefVect/norm(RefVect); 

CosThetaObjPrime(i,j)=  dot(ObjNormalHat,RefV  ectHat); 
end 
end 


% . - . - . — . -- . — . - 

%Cos  of  the  angle  between  the  Imaged  Wall  Normal  and  every  vector  between 
%the  Object  Points  and  thelmaged  Wall  points 

CosThetaWalllm  =  zeros(TotObjPoints,TotWallImPoints); 
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for  i=l:TotObjPoints 
for  j=l  :TotWallImPoints 

IncVect  =  FlatObjArray(i,:)-FlatWallImArray(j,:); 
IncVectHat  =  IncVect/norm(IncVect); 
CosThetaWallIm(i,j)=  dot(WallNormalHat, IncVectHat); 
end 
end 


% . — . - . - . 

%Cos  of  the  angle  between  the  Wall  Normal  and  every  vector  between  the  Wall 
%points  and  the  Lens  Points 

CosThetaWalllmPrime  =  zeros(TotWallImPoints,TotLensPoints); 

for  i=l:TotWallImPomts 
for  j=l:TotLensPoints 

RefVect  =  FlatLensArray(j,:)-FlatWallImArray(i,:); 

RefVectHat  =  RefVect/norm(RefVect); 

CosThetaWallImPrime(i,j)  =  dot(WallNormalHat, RefVectHat); 
end 
end 


% - 

%Cos  of  the  angle  between  the  Lens  Normal  and  every  vector  between  the  Wall 
%points  and  the  Lens  Points 

CosThetaLens  =  zeros(TotWallImPoints,TotLensPoints); 

for  i=l:TotWallImPoints 
for  j=l  :TotLensPoints 

IncVect  =  FlatWallImArray(i,:)-FlatLensArray(j,:); 

IncVectHat  =  IncVect/norm(IncVect); 

CosThetaLens(i,j  )=  dot(LensNormalHat,Inc  V  ectHat) ; 
end 
end 

clear  (’ Wall  Array  VWalllmArray’) 


% . . - . — . — . - .  . 

%  This  section  creates  the  Glint  Anlge  Arrays  for  calculations 
% - 


% - - - - - 

%  This  section  creates  the  Wall  Glint  Angle  Array 

WallGlintAngle=  zeros(TotWallPoints,TotObjPoints); 

IncRay  =  laser; 
for  i  =  1  :TotWallPoints 
for  j  =l:TotObjPoints 

RefRay  =F  latObj  Array(j , : )  -FlatW  all  Array(i, : ) ; 
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RefRayHat  =  RefRay/norm(RefRay); 

GlintVec  =  (RefRayHat  +  IncRay)/2; 

GlintVecHat  =  Glint  Vec/norm(GlintVec); 
WallGhntAngle(ij)=acos(dot(GlintVecHat,WallNormalHat)); 
end 
end 


%  This  section  creates  the  Object  Glint  Angle  Array 

ObjGlintAngle  =  zeros(TotWallPoints,TotObjPoints,TotWallImPoints); 

for  i  =  1  :TotWallPoints 
for  j=l:TotObjPoints 

IncRay  =  FlatWallArray(i,:)-FlatObjArray(j,:); 

IncRayHat  =  IncRay/norm(IncRay); 
for  k  =  1  :TotWallImPoints 

RefRay  =  FlatWallImArray(k,:)-FlatObjArray(j,:); 

RefRayHat  =  RefRay/norm(RefRay); 

GlintVec  =  (RefRayHat  +  IncRayHat)/2; 

GlintVecHat  =  GlintVec/norm(GlintVec); 

Obj  Glint  Angle(i,j,k)=  acos(dot(GlintVecHat,ObjNormalHat)); 
end 
end 
end 


% - - - — . — . ----- 

%  This  section  creates  the  Wall  Imaged  Glint  Angle  Array 

WalllmGlint  Angle  =  zeros(TotObjPoints,TotWallImPoints,TotLensPoints); 

for  i  =  1  :TotObj  Points 
for  j  =  HTotWalllmPoints 

IncRay  =  FlatObjArray(i,:)-FlatWallImArray(j,:); 

IncRayHat  =  IncRay/norm(IncRay); 
for  k  =  1  :TotLensPoints 

RefRay  =  FlatLensArray(k,:)-FlatWallImArray(j,:); 

RefRayHat  =  RefRay/norm(RefRay); 

GlintVec  =  (RefRayHat  +  IncRayHat)/2; 

GlintVecHat  =  GlintVec/norm(GlintVec); 

WalllmGlint  Angle(i,j  ,k)=  acos(dot(GlintVecHat,WallNormalHat)); 
end 
end 
end 

clear  (’FlatLensArrayVFlatObjArrayVFlatWallArrayVFlatWalllmArray’) 
clear  (’RefRay VRefRayHaf,’RefV ecf,’RefV ectHaf ) 
clear  (’ObjNormalHaf,’LensNormalHatVWallNormalHaf) 
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% - - - - 

%  This  section  creates  the  model 
% . . - - - 


% . . . . 

%  This  section  models  the  Wall  Radiance 

LWall  =  zeros(TotWallPoints,TotObjPoints); 
for  i  =  1  :TotWallPoints 
for  j  =  l:TotObjPoints 

LWall(i,j)  =(l/(2*pi*SigmaWallA2*cos(WallGlintAngle(i,j))A3))*... 
exp(-tan(WallGlintAngle(i,j)A2)/(2*  SigmaWallA2)); 
end 
end 


% - 

%  This  section  models  the  Object  Irradiance 

EObj  =  zeros(TotWallPoints,TotObjPoints); 

for  i  =  1  :TotWallPoints 
for  j = 1 :  T  otObj  P  oints 
EObj(i,j)=  LWall(ij)*... 

CosThetaWallPrime(i,j)*CosThetaObj(i,j)/DistToObj(i,j)A2; 

end 

end 

WGA  =  reshape(WallGlintAngle(5,:),xNumObjPoints,zNumObjPoints); 
LW  =  reshape(LWall(5,:),xNumObjPoints,zNumObjPoints); 

EO  =  reshape(LWall(5,:),xNumObjPoints,zNumObjPoints); 

clear  (’WallGlintAngleVCosThetaWallPrimeVCosThetaObj’,... 

’DistT  oObj  ’,’LW  all’) 


% - 

%  This  section  models  the  Object  Radiance 

LObj  =  zeros(TotWallPoints,TotObjPoints,TotWallImPoints); 

for  i  =  1  :TotWallPoints 
for  j  =  ETotObjPoints 
for  k  =  1  :TotWallImPoints 
LObj(i,j,k)=  EObj(i,j)*... 

( 1  /(2  *pi  *  SigmaObj  A2  *cos(Obj  Glint  Angle(i  j  ,k))  A3 ))  * . . . 
exp(-tan(ObjGlintAngle(i,j,k)A2)/(2*  SigmaObj A2)); 
end 
end 
end 
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% - - - — . — - - -  - 

%  This  section  models  the  Imaged  Wall  Irradiance 

Elm  =  zeros(TotWallPoints,TotObjPoints,TotWallImPoints); 

for  i  =  1  :TotWallPoints 
for  j  =  ETotObjPoints 
for  k  =  1  :TotWallImPoints 
EIm(i,j,k)=  LObj(i,j,k)*... 

CosThetaObjPrimeO,k)*CosThetaWallImO’,k)/DistToWallImO,k)A2; 

end 

end 

end 

clear(’EObj ’,’Obj  GlintAngleVLObj ’,’CosThetaObj  Prime’, .. . 

’CosThetaW  alllm’,' DistT  oW  alllm’) 


% .  . . . - .  . . 

%  This  section  models  the  Imaged  Wall  Radiance 

LIm  =zeros(TotWallPoints,TotObjPoints,TotWallImPoints, TotLensPoints); 

for  i  =  1  :TotWallPoints 
for  j  =  ETotObjPoints 
for  k  =  1  :TotWallImPoints 
for  1  =  1  :TotLensPoints 
LIm(i,j,k,l)=  EIm(i,j,k)*... 

(l/(2*pi*SigmaWallA2*cos(WallImGlintAngle(j,k,l))A3))*... 
exp(-tan(WallImGlintAngle(j,k,l)A2)/(2*  SigmaWallA2)); 
end 
end 
end 
end 


% . . — - - - — . — . . 

%  This  section  models  the  Imaged  Wall  Radiance 

ELens=  zeros(TotWallPoints,TotObjPoints,TotWallImPoints,... 
TotLensPoints); 

for  i  =  1  :TotWallPoints 
for  j  =  ETotObjPoints 
for  k  =  1  :TotWallImPoints 
for  1  =  1:  TotLensPoints 

ELens(i,j,k,l)=  LIm(i,j,k,l)*  ... 

CosThetaWallImPrime(k,l)  *CosThetaLens(k,l)/. . . 
DistToLens(k,l)  A2 ; 
end 
end 
end 
end 
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clear  (’WalllmGlint Angle VEImVCosThetaWalllmPrimeVCosThetaLens’,... 
’DistT  oLens ’,’LIm’) 

clear  ('GlintVecVGlintVecHat'/IncRayVIncRayHatVIncVect’,... 
’IncVectHat’) 


% - 

%  This  section  creates  a  2D  dual  image  of  the  center  point  object 

DualImage2D  =  zeros(21,21); 

for  x  =  3:23 

xSlideEdgePlus  =  x+2; 
xSlideEdgeMinus  =  x-2; 

for  z  =  3:23 

zSlideEdgePlus  =  z+2; 
zSlideEdgeMinus  =  z-2; 

Slide  =  zeros(xNumObj  Points, zNumObj  Points); 

Slide(x,z)=l; 

FlatSlide  =  squeeze(reshape(Slide,l,TotObjPoints)); 

for  i  =  ETotWallPoints 
ilium  =  1 ; 
switch  i 
case  2 
ilium  =  3; 
case  4 
ilium  =  3; 
case  5 
ilium  =  5; 
case  6 
ilium  =  3; 
case  8 
ilium  =  3; 
end 

for  jl  =  11:13 
for  j2  =  12:14 
j  =  25*jl+j2; 

for  k  =  1  :TotWallImPoints 

for  1  =  (TotLensPoints+l)/2:(TotLensPoints+l)/2 
temp  =  ELens(i,j,k,l)*FlatSlide(j); 
DualImage2D(x-2,z-2)=. . . 
DualImage2D(x-2,z-2)+temp; 
end 
end 
end 
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end 

end 

end 

end 

DPoint  =  DualImage2D; 


% - 

%  This  section  creates  a  2D  indirect  image  of  the  center  point  object 

IndirectImage2D  =  zeros(21,21); 

for  x  =  3:23 

xSlideEdgePlus  =  x+2; 
xSlideEdgeMinus  =  x-2; 

for  z  =  3:23 

zSlideEdgePlus  =  z+2; 
zSlideEdgeMinus  =  z-2; 

Slide  =  zeros(xNumObjPoints,zNumObjPoints); 

Slide(x,z)=l; 

FlatSlide  =  squeeze(reshape(Slide,l,TotObjPoints)); 

for  i  =  ETotWallPoints 
ilium  =  1 ; 
switch  i 
case  2 
ilium  =  3; 
case  4 
ilium  =  3; 
case  5 
ilium  =  5; 
case  6 
ilium  =  3; 
case  8 
ilium  =  3; 
end 

for  j  =  ETotObjPoints 
for  k  =  1  :TotWallImPoints 

for  1  =  (TotLensPoints+l)/2:(TotLensPoints+l)/2 
temp  =  ELens(i,j,k,l)*FlatSlide(j); 
IndirectImage2D(x-2,z-2)=. . . 
IndirectImage2D(x-2,z-2)+temp; 
end 
end 
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end 


end 

end 

end 


IPoint  =  IndirectImage2D; 


% . - — . — . — . . . . 

%  This  section  creates  a  2D  dual  image  of  the  two  squares 

DualImage2D  =  zeros(21,21); 

for  x  =  3:23 

xSlideEdgePlus  =  x+2; 
xSlideEdgeMinus  =  x-2; 

for  z  =  3:23 

zSlideEdgePlus  =  z+2; 
zSlideEdgeMinus  =  z-2; 

Slide  =  zeros(xNumObjPoints,zNumObjPoints); 

Slide(x-l,z-l)=l; 

Slide(x+l,z+l)=l; 

FlatSlide  =  squeeze(reshape(Slide,l,TotObjPoints)); 

for  i  =  ETotWallPoints 
ilium  =  1 ; 
switch  i 
case  2 
ilium  =  3; 
case  4 
ilium  =  3; 
case  5 
ilium  =  5; 
case  6 
ilium  =  3; 
case  8 
ilium  =  3; 
end 

for  jl  =  11:13 
forj2=  12:14 
j=25*jl+j2; 
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for  k  =  1  :TotWallImPoints 

for  1  =  (TotLensPoints+l)/2:(TotLensPoints+l)/2 
temp  =  ELens(i,j,k,l)*FlatSlide(j); 
DualImage2D(x-2,z-2)=. . . 
DualImage2D(x-2,z-2)+temp; 
end 
end 
end 

end 

end 

end 

end 

D2Point  =  DualImage2D; 


% .  . - . . . . 

%  This  section  creates  a  2D  indirect  image  of  the  two  squares 

IndirectImage2D  =  zeros(21,21); 

for  x  =  3:23 

xSlideEdgePlus  =  x+2; 
xSlideEdgeMinus  =  x-2; 

for  z  =  3:23 

zSlideEdgePlus  =  z+2; 
zSlideEdgeMinus  =  z-2; 

Slide  =  zeros(xNumObjPoints,zNumObjPoints); 

Slide(x-l,z-l)=l; 

Slide(x+l,z+l)=l; 

FlatSlide  =  squeeze(reshape(Slide,l,TotObjPoints)); 

for  i  =  ETotWallPoints 
ilium  =  1 ; 
switch  i 
case  2 
ilium  =  3; 
case  4 
ilium  =  3; 
case  5 
ilium  =  5; 
case  6 
ilium  =  3; 
case  8 
ilium  =  3; 
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end 

for  j  =  l:TotObj  Points 
for  k  =  1  :TotWallImPoints 

for  1  =  (TotLensPoints+l)/2:(TotLensPoints+l)/2 
temp  =  ELens(i,j,k,l)*FlatSlide(j); 
IndirectImage2D(x-2,z-2)=. . . 
IndirectImage2D(x-2,z-2)+temp; 
end 
end 
end 


end 


end 

end 


I2Point  =  IndirectImage2D; 
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Appendix  C.  Reflector  BRDF  measurements 


Figure  62  shows  the  measured  BRDF  from  the  semi-gloss  reflector  resulting  from 
633nm  HeNe  laser  and  a  45  degree  incident  angle.  In  this  graph,  the  ‘X’  axis  is  the  angle 
difference  from  the  specular  reflection  i.e.  0  on  the  ‘X’  axis  represents  specular  reflect, 
positive  angles  are  away  from  the  surface  normal  and  negative  angles  are  back  toward  the 
surface  normal  and  the  incident  irradiance.  The  gap  in  the  data  at  -90  degrees  comes 
from  the  sensor  blocking  the  incident  irradiance.  The  data  was  gathered  from  three 
different  locations  to  validate  the  assumption  the  reflector  was  homogenous.  In  one  of 
the  measurements,  the  sample  was  rotated  90  degrees,  to  validate  the  isotropic 
assumption. 
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Figure  62.  Measured  semi-gloss  reflector  BRDF 


Figure  63  shows  the  measured  BRDF  from  the  semi-gloss  reflector  resulting  from 
633nm  HeNe  laser  and  a  45  degree  incident  angle. 
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Figure  63.  Measured  flat  white  reflector  BRDF 
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Appendix  D.  Semi-gloss  reflector  1-D  images 


Figure  32,  repeated  here  as  Figure  64,  shows  the  MTF  created  from  the  1-D 
experiment  using  the  semi-gloss  imaged  reflector.  The  raw  and  deconvolved  dual  and 
indirect  images  as  well  as  the  Fourier  transforms  used  to  create  the  MTF  are  shown  in 
Figure  65  through  69. 
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Figure  64.  1-D  Semi-gloss  MTF 


Figure  65  shows  the  dual  images  created  using  the  semi-gloss  imaged  reflector  of 
the  following  slides:  (a)  0.1  cycles/mm,  (b)  0.2  cycles/mm,  (c)  0.5  cycles./mm,  (d)  1.0 
cycles/mm,  (e)  2.0  cycles/mm  and  (f)  3.0  cycle/mm.  Figure  65  (g),  (h)  and  (i)  show  an 
expanded  view  of  (d),  (e)  and  (f)  respectively.  Figure  66  shows  the  1-D  indirect  images 
created  using  the  semi-gloss  imaged  reflector  of  the  following  slides:  (a)  0. 1  cycles/mm, 
(b)  0.2  cycles/mm,  (c)  0.5  cycles/mm,  (d)  1.0  cycles/mm,  (e)  2.0  cycles/mm  and  (f)  3.0 
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cycle/mm.  Figure  66  (g),  (h)  and  (i)  show  an  expanded  view  of  (d),  (e)  and  (f) 
respectively.  Figure  67  shows  the  deconvolved  images  overlaid  on  the  raw  images,  (a) 
through  (f)  are  the  dual  images  while  (g)  through  (1)  are  the  indirect  images.  The 
individual  images  are  as  follows:  (a)  and  (g)  0.1  cycles/mm,  (b)  and  (h)  0.2  cycles/mm, 
(c)  and  (i)  0.5  cycles/mm,  (d)  and  (j)  expanded  view  of  1.0  cycles/mm,  (e)  and  (k) 
expanded  view  of  2.0  cycles/mm,  (f)  and  (1)  expanded  view  of  the  3.0  cycle/mm.  Figure 
68  (a)  through  (1)  shows  the  Fourier  transforms  of  the  dual  images  shown  in  Figure  67  (a) 
through  (1),  while  Figure  69(a)  through  (1)  shows  the  Fourier  transforms  of  the  indirect 
images  of  Figure  67  (a)  through  (1). 


116 


Mean  Intensity  Mean  Intensity  Mean  Intensity 


2000 


^  500 


°0  10  20  30  40  50  60  70  80  90  100 

milli-meters 


(a) 


(b) 


(c) 


0  10  20  30  40  50  60  70  80  90  100 

milli-meters 


milli-meters 


(d) 


(e) 


(f) 


47  49  51  53 

milli-meters 


(g) 


(h) 


(0 


Figure  65.  1-D  Semi-gloss  dual  images 
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Figure  66.  1-D  Semi-gloss  indirect  images 
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Figure  67.  1-D  Semi-gloss  dual  and  indirect  deconvolved  images 
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Figure  68.  1-D  Semi-gloss  dual  and  indirect  image  Fourier  transforms 
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Figure  69.  1-D  Semi-gloss 
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deconvolved  dual  and  indirect  image  Fourier  transforms 
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Appendix  E.  Flat  White  reflector  1-D  images 


Figure  33,  repeated  here  as  Figure  70,  shows  the  MTF  created  from  the  1-D 
experiment  using  the  flat  white  imaged  reflector.  The  raw  and  deconvolved  dual  and 
indirect  images  as  well  as  the  Fourier  transforms  used  to  create  the  MTF  are  shown  in 
Figure  71  through  75. 
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Figure  70.  1-D  Flat  white  MTF 

Figure  71  shows  the  dual  images  created  using  the  flat  white  imaged  reflector  of 
the  following  slides:  (a)  0.1  cycles/mm,  (b)  0.2  cycles/mm,  (c)  0.5  cycles./mm,  (d)  1.0 
cycles/mm,  (e)  2.0  cycles/mm  and  (f)  3.0  cycle/mm.  Figure  71  (g),  (h)  and  (i)  show  an 
expanded  view  of  (d),  (e)  and  (f)  respectively.  Figure  72  shows  the  1-D  indirect  images 
created  using  the  flat  white  imaged  reflector  of  the  following  slides:  (a)  0.1  cycles/mm, 
(b)  0.2  cycles/mm,  (c)  0.5  cycles/mm,  (d)  1.0  cycles/mm,  (e)  2.0  cycles/mm  and  (f)  3.0 
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cycle/mm.  Figure  72  (g),  (h)  and  (i)  show  an  expanded  view  of  (d),  (e)  and  (f) 
respectively.  Figure  73  shows  the  deconvolved  images  overlaid  on  the  raw  images  (a) 
through  (f)  are  the  dual  images  while  (g)  through  (1)  are  the  indirect  images.  The 
individual  images  are  as  follows:  (a)  and  (g)  0.1  cycles/mm,  (b)  and  (h)  0.2  cycles/mm, 
(c)  and  (i)  0.5  cycles/mm,  (d)  and  (j)  expanded  view  of  1.0  cycles/mm,  (e)  and  (k) 
expanded  view  of  2.0  cycles/mm,  (f)  and  (1)  expanded  view  of  the  3.0  cycle/mm.  Figure 
74  (a)  through  (1)  shows  the  Fourier  transforms  of  the  dual  images  shown  in  Figure  73  (a) 
through  (1),  while  Figure  75(a)  through  (1)  shows  the  Fourier  transforms  of  the  indirect 
images  of  Figure  73  (a)  through  (1). 
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Figure  71.  1-D  Flat  white  dual  images 
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Figure  72.  1-D  Flat  white  indirect  images 
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Figure  73.  1-D  Flat  white  dual  and  indirect  deconvolved  images 
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Figure  74.  1-D  Flat  white  dual  and  indirect  image  Fourier  transforms 
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Figure  75.  1-D  Flat  white  deconvolved  dual  and  indirect  image  Fourier  transforms 
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Appendix  F.  Unpainted  reflector  1-D  images 


Figure  34,  repeated  here  as  Figure  76,  shows  the  MTF  created  from  1-D 


experiment  using  the  unpainted  imaged  reflector.  The  raw  and  deconvolved  dual  and 


indirect  images  as  well  as  the  Fourier  transforms  used  to  create  the  MTF  are  shown  in 


Figure  77  through  81. 
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Figure  76.  1-D  Unpainted  MTF 


Figure  77  shows  the  dual  images  created  using  the  unpainted  imaged  reflector  of 


the  following  slides:  (a)  0.1  cycles/mm,  (b)  0.2  cycles/mm,  (c)  0.5  cycles./mm,  (d)  1.0 


cycles/mm,  (e)  2.0  cycles/mm  and  (f)  3.0  cycle/mm.  Figure  65  (g),  (h)  and  (i)  show  an 


expanded  view  of  (d),  (e)  and  (f)  respectively.  Figure  78  shows  the  1-D  indirect  images 


created  using  the  unpainted  imaged  reflector  of  the  following  slides:  (a)  0.1  cycles/mm, 


(b)  0.2  cycles/mm,  (c)  0.5  cycles/mm,  (d)  1.0  cycles/mm,  (e)  2.0  cycles/mm  and  (f)  3.0 
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cycle/mm.  Figure  78  (g),  (h)  and  (i)  show  an  expanded  view  of  (d),  (e)  and  (f) 
respectively.  Figure  79  shows  the  deconvolved  images  overlaid  on  the  raw  images,  (a) 
through  (f)  are  the  dual  images  while  (g)  through  (1)  are  the  indirect  images.  The 
individual  images  are  as  follows:  (a)  and  (g)  0.1  cycles/mm,  (b)  and  (h)  0.2  cycles/mm, 
(c)  and  (i)  0.5  cycles/mm,  (d)  and  (j)  expanded  view  of  1.0  cycles/mm,  (e)  and  (k) 
expanded  view  of  2.0  cycles/mm,  (f)  and  (1)  expanded  view  of  the  3.0  cycle/mm.  Figure 
80  (a)  through  (1)  shows  the  Fourier  transforms  of  the  dual  images  shown  in  Figure  79  (a) 
through  (1),  while  Figure  81(a)  through  (1)  shows  the  Fourier  transforms  of  the  indirect 
images  of  Figure  79  (a)  through  (1). 
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Figure  77.  1-D  Unpainted  dual  images 
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Figure  78.  1-D  Unpainted  indirect  images 
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Figure  79.  1-D  Unpainted  dual  and  indirect  deconvolved  images 
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Figure  80.  1-D  Unpainted  dual  and  indirect  image  Fourier  transforms 
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Figure  81.  1-D  Unpainted  deconvolved  dual  and  indirect 
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